UNCLASSIFIED 


AD  NUMBER 


AD485025 


NEW  LIMITATION  CHANGE 
TO 

Approved  for  public  release,  distribution 
unlimited 


FROM 

Distribution  authorized  to  U.S.  Gov't, 
agencies  and  their  contractors; 
Administrative/Operational  Use;  JUL  1966. 
Other  requests  shall  be  referred  to  Air 
Force  Weapons  Lab.,  Kirtland  AFB,  NM. 


AUTHORITY 


AFWL  per  DTIC  form  55 


THIS  PAGE  IS  UNCLASSIFIED 


485025 


AFWL-TR-65-143,  Vol  IV 


AFWL-TR 
65-143 
Vol  IV 


FIREBALL  PHENOMENOLOGY 
AND  CODE  DEVELOPMENT 


Volume  IV 


SPUTTER  Subroutines  for  Radiation  Transport  in  Planes 


General  Atomic  Division  of  General  Dynamics  Corporation 
Special  Nuclear  Effects  Laboratory 
San  Diego,  California 
Contract  AF  29(601)-6492 


TECHNICAL  REPORT  NO.  AFWL-TR-65-143,  Vol  IV 


July  1966 


V 


AIR  FORCE  WEAPONS  LABORATORY 
Research  and  Technology  Division 
Air  Force  Systems  Command 
Kir'land  Air  Force  Base 
New  Mexico 


AFWL-'-  F.  *65-143,  Vol  IV 


Research  and  Technology  Division 
AIR  FORCE  WEAPONS  LABORATORY 
Air  Force  Systems  Command 
Kirtland  Air  Force  Base 
New  Mexico 


When  U.  S.  Government  drawings,  specifications,  or  other  data  are  used  for 
any  purpose  other  than  a  definitely  related  Government  procurement  operation, 
the  Government  thereby  incurs  no  responsibility  nor  any  obligation  whatsoever, 
and  the  fact  that  the  Government  may  have  formulated,  furnished,  or  in  any 
way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not  to  be 
regarded  by  implication  or  otherwise,  as  in  any  manner  licensing  the  holder 
or  any  other  person  or  corporation,  or  conveying  any  rights  or  permission  to 
manufacture,  use,  or  sell  any  patented  invention  that  may  in  any  way  be 
related  thereto. 

This  report  is  made  available  for  study  with  the  understanding  that 
proprietary  interests  in  and  relating  thereto  will  not  be  impaired.  In 
case  of  apparent  conflict  or  any  other  questions  between  the  Government's 
rights  and  those  of  others,  notify  the  Judge  Advocate,  Air  Force  Systems 
Command,  Andrews  Air  Force  Base,  Washington,  D.  C.  20331. 

This  document  is  subject  to  special  export  controls  and  each  transmittal 
to  foreign  governments  or  foreign  nationals  may  be  made  only  with  prior 
approval  of  AFWL  (WLRT) ,  Kirtland  AFB,  NM,  87117.  Distribution  is  limited 
because  of  the  technology  discussed  in  the  report. 


/ 


AFWL-TR-65-143,  Vol  IV 


FIREBALL  PHENOMENOLOGY  AND  CODE  DEVELOPMENT 
Volume  IV 

SPUTTER  Subroutines  for  Radiation  Transport  in  Planes 

General  Atomic  Division  o£  General  Dynamics  Corporation 
Special  Nuclear  Effects  Laboratory 
San  Diego,  California 
Conti act  AF  29 (601) -6492 


TECHNICAL  REPORT  NO.  AFWL-TR-65-143,  Vol  IV 


This  document  is  subject  to  special 
export  controls  and  each  transmittal  to 
foreign  governments  or  foreign  nationals 
may  be  made  only  with  prior  approval  of 
AFWL  (WLRT),  Kirtland  AFB,  NM,  87117. 
Di8tributiou  is  limited  because  of  the 
technology  discussed  in  the  report. 


FOREWORD 


AFWL-TR-65-143  ^  Vol  IV 


This  report  was  prepared  by  General  Atomic  Division  of  General  Dynamics 
Corporation,  San  Diego,  California,  under  Contract  AP  29 (601) -6492.  The 
research  was  performed  under  Program  Element  7. 60, 06. 01, Df,  Project  5710, 

Subtask  07.003/005,  and  was  funded  by  the  Defense  Atomic  Support  Agency  (DASA) . 

Inclusive  dates  of  research  were  1  June  1963  to  13  July  1965.  The  report 
was  submitted  15  March  1966  by  the  Air  Force  Weapons  Laboratory  Project  Officer, 
lLt  F.  C.  Tompkins  III  (WLRT) .  The  contractor’s,  report  number  is  GA-6585. 

This  report  is  divided  ir.to  six  volumes  as  follows:  Volume  I,  Summary  and 
the  Fireball  Models;  Volume  II,  Early  Fireball  Phenomena  in  the  TIGHTROPE  Event;* 
Volume  III,  SPUTTER  Subroutines  for  Radiation  Transport  in  Spheres*;  Volume  IV, 
SPUTTER  Subroutines  for  Radiation  Transport  in  Planes*;  Volume  V,  Material 
Properties;  and  Volume  VI,  Extensions  of  the  Physics  and  Problem  Areas. 

The  SPUTTER  subroutines  for  radiation  transport  in  planes  described  in 
Volume  IV  were  developed  by  Dr.  B.  E.  Freeman  and  Dr.  C.  G.  Davis,  Jr.  The 
cooperation  and  contributions  of  Captains  Milton  Gillespie,  William  Whittaker, 
and  George  Spillman  of  the  Air  Force  Weapons  Laboratory  are  gratefully 
acknowledged. 

This  technical  report  has  been  reviewed  and  is  approved. 


Project  Officer 


RALPH  /.  PENT.INGTON  f 

Lt  Colonel,  USAF 

Chief,  Theoretical  Branch 


Chief,  Research  Division 


*  Volume  II  has  been  withdrawn  ,  and  will  not  be  published. 


ii 


AFWL-TR-65-143,  Vol  IV 


ABSTRACT 


The  radiation-transport  subroutines  of  the  SPUTTER  code  for  plane  slab  geometry 
have  been  supplemented  by  alternative,  formulation  based  on  integration  along 
sampling  rry  paths  through  the  slab.  Angular  Integrations  are  performed  by 
the  Gaussian  quadrature  method  which  determines  the  ray  angles.  Options  may 
be  exercised  to  determine  the  number  of  angles  and  the  nature  of  the  radiation 
boundary  condition  at  one  boundary  of  the  transport  region.  The  characteristic 
ray  code  differ i  from  the  current  integral  method  in  performing  problems  having 
a  large  number  of  zones  more  rapidly  and  in  having  more  general  boundary 
conditions.  For  most  applications  a  small  number  of  angles  give  adequate 
accuracy.  The  numerical  method  used  in  the  ray  code  is  described.  In  addition, 
the  organization  of  the  code  is  discussed  end  subroutines  are  listed. 
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* 


The  SPUTTER  code  subroutines  for  radiation  transport  in  planes  described 
herein  are  as  they  existed  on  July  30,  1965.  General  Atomic  has  exercised  due 
care  in  preparation,  but  does  not  warrant  the  merchantability,  accuracy,  and 
completeness  of  these  subroutines  or  of  their  description  contained  herein. 

The  complexity  of  this  kind  of  program  precludes  any  guarantee  to  that  effect. 
Therefore,  any  user  must  make  his  own  determination  of  the  suitability  of  these 
subroutines  for  any  specific  \ise  and  of  the  validity  of  the  information  produced 
by  their  use. 
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SECTION  1 


INTRODUCTION 


The  new  routines  for  radiation  transport  in  planes  closely  parallel  the 
spherical  radiation  transport  subroutines  (Volume  III)  in  mathematical 
formulation  and  code  organization.  This  parallelism  is  especially  close 
between  the  SRADTN  (for  spheres)  and  PRADTN  (for  planes)  subroutines  in 
which  calculations  peripheral  to  the  transport  integration  are  performed. 

.  In  fact,  it  is  likely  that  these  subroutines  can  be  condensed  to  form  a  single 
subroutine  for  spheres  and  planes,  although  currently  they  are  separate. 

The  routines  reported  here  are  to  be  considered  as  alternatives  to 
the  routines  based  on  the  integral  formulation  of  the  transport  equation 
currently  n  use.  By  comparison,  the  current  integral  version  is  more 
accurate  in  the  performance  of  angular  integrations  of  the  intensity,  but 
for  problems  requiring  a  large  number  of  zones  it  requires  more  com¬ 
putation  time.  Boundary  conditions  on  the  radiation  intensity,  however, 
are  much  more  naturally  incorporated  into  the  new  version. 

Conditions  suggesting  a  preference  for  using  the  new  routines  are: 

(1)  a  large  number  of  zones  and  a  desire  to  reduce  calculation  time  and 

(2)  the  necessity  of  specifying  a  radiation  intensity  incident  on  the  slab 
surface  which  has  arbitrary  angular  and  frequency  dependence. 

The  numerical  sequences  used  in  solving  the  transport  equation  are 
discussed  in  Section  II.  A  description  of  the  diffusion  approximation  used 
in  conjunction  with  the  transport  solution  is  given  in  Section  III  and  a  brief 
description  of  the  methods  of  frequency  integration  is  given  in  Section  IV. 
Section  V  includes  the  actual  code  description  in  terms  of  code  organiza¬ 
tion  and  economics.  Section  VI  includes  some  initial  studies  on  timing 
and  accuracy  in  the  angular  integrations. 
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SECTION  II 

NUMERICAL  SOLUTION  OF  THE  TRANSPORT  EQUATION 


The  radiation  routines  described  herein  contain  a  formulation  based 
on  the  numerical  solution  of  the  radiation  transport  equation  along  a  selec¬ 
tion  of  sampling  rays  through  the  slab.  Relevant  averages  over  the 
angular  distribution  are  obtained  by  numerical  quadrature,  as  described 
in  Section  2.  3,  and  the  numerical  solution  of  the  transport  equation  along 
the  photon  ray  is  presented  in  Section  2.  1.  Criteria  for  selecting  the 
sampling  rays  are  discussed  in  Section  2.  2.  All  of  the  derivations  of  this 
section  apply  to  photons  of  a  particular  frequency;  integration  over  fre¬ 
quency  is  discussed  in  Section  IV. 


The  radiation  transport  equation  in  plane  geometry  that  describes 
the  changes  in  the  .specific  intensity  Iy  of  photons  of  frequency  v  resulting 
from  pure  absorption  and  emission  according  to  the  local  thermodynamic 
equilibrium  assumption  is 
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(2.1) 


where 
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and  <rv  is  the  pure  absorption  coefficient.  The  scattering  coefficient  is 
assumed  to  be  negligibly  small  compared  to  the  absorption  coefficient. 
Additionally,  the  retardation  of  the  photons  is  neglected,  as  is  valid  when 
the  radia^on  energy  is  small  and  temperatures  change  slowly.  The 
resulting  equation  descrioes  the  qvasi- steady  intensity  field  resulting  from 
the  distribution  of  sources  existing  at  a  particular  time. 


Defining  the  monochromatic  optical  depth,  t,  as 
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For  the  case  of  a  constant  or  step- function  source,  the  source  function  B 
takes  a  value  dependent  on  which  interface  of  the  zone  is  affected.  If  the 
left  interface  (r  =  r.  satisfies  the  criteria  for  a  constant  source, 

B  =  B.  i,  t.  ,  <‘t  <  r.  i  . 
i-2  i-l  i-i 

If  the  right  interface  (t  =  r^)  satisfies  the  criteria, 

i-£  i-f  “  1 


The  integral  of  Eq.  (2.  3)  can  be  evaluated  with  the  interpolation  function 
of  Eq.  (2.  4)  to  give  for  the  intensity 
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where 
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In  these  expressions,  A  =  -  r^.  The  coefficients  of  Eq.  (2,  5)  can  be 

re-expressed  by  using  the  definitions  of  Eq.  (2.  4): 
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The  terms  in  Eq.  (2.  6)  may  be  interpreted  as  containing  combinations  of 
numerical  approximations  to  the  values  o'f  the  source  function  and  the  t 
derivative  of  the  source  function  at  the  boundaries  of  the  interval. 


This  form  of  the  equation,  in  fact,  can  be  obtained  in  another  way 
starting  from  Eq.  (2.  3).  Two  successive  integrations  by  parts  transforms 
the  expression  for  I.  into  the  following  equivalent  form: 


-(T 


1 


-T) 


e 


dr 


9 


(2.7) 


in  terms  of  values  of  the  source  function  and  the  first  two  derivatives  of 
the-source  function  with  respect  to  t. 


In  an  optically  thin  interval,  the  most  important  contribution  arises 
from  the  termo  I^_i  and  B,  which  represent  the  transmitted  intensity  and 
the  emission  from  the  zone.  The  derivative  terms  cancel  in  this  approx¬ 
imation;  this  is  perhaps  more  directly  indicated  by  Eq.  (2.  3).  In  the 
optically  thick  interval,  which  is  the  extreme  opposite,  only  the  first  two 
terms  evaluated  at  i  are  usually  of  significance.  The  terms  £rom  i  -  1 
are  strongly  attenur  ;ed  and  8^B/9t^  in  the  integral  is  usually  small.  In 
the  limit,  the  diffusion  approximation  results  from  the  term  9B/9T)i. 
Between  limits,  it  is  necessary  to  consider  the  integral  term  in  Eq.  (2.  7). 


If  A  is  not  too  large,  a  representative  mean  value  of  the  exponential 
in  the  interval  may  be  taken  to  give  for  the  integral  of  Eq.  (2.  7) 
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and  thus  the  expression  for  intensity  becomes 
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This  expression  has  just  the  form  of  Eqs.  (2.  5)  and  (2.  6)  when  the  dif¬ 
ference  expressions  are  identified  with  the  derivatives, 


It  is  clear  from  the  derivation  of  Eq.  (2.  5)  that  the  resulting  intensity 
is  a  positive  quantity.  With  positive  values  for  zone  source  functions,  the 
linear  interpolation  expression  assures  that  the  integral  contribution  is 
always  positive.  Since  the  boundary  intensity  is  alway.3  a  positive  quantity, 
the  positivity  of  all  intensities  is  assured. 


In  the  diffusion  approximation  limit,  only  quantities  at  interface  i  will 
survive,  and 


which  can  be  evaluated  as 
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The  independent  variable  h  depends  only  on  x,  so  that  angular  integrations 
of  I^  can  be  performed  explicitly  in  the  diffusion  approximation,  which 
takes  account  of  the  dependence  on  angle  of  Eq.  (2.  9).  A  difference 
approximation  can  also  be  based  on  this  expression,  assuming  that  B  is 
linear  in  h,  i.  e. , 
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where  fi  is  the  cosine  of  the  angle  which  the  ray  makes  with  the  slab 
normal.  The  corresponding  equation  for  the  intensity  is  Eq.  (2.  5),  in 
which 
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2.  1.  2.  Small-optical- depth  Expansion  , 

If  the  optical  depth  is  very  small,  the  intensity. expression.tin 
Eq.  (2.  3)  takes  a  much  simpler  form, 
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Although  this  result  is  the  limiting  form  of  Eqs,  (2,  5)  and  (2.  6),  the  terms 
must  cancel  through  second  order  in  an  expansion  in.  A  before  the  first 
surviving  term,  derived  in  part  from  the  quadratic  terms  of  the  exponen¬ 
tials,  is  obtained.  Consequently,  for  sufficiently  small  argument,  the 
finite  number  of  figures  used  in  the  exponential  will  render  the  result 
inaccurate.  For  the  exponential  from  the  IBM-7044  system,  this  restricts 
che  argument  to  a  number  greater  than  ~  2  x  10-4;  but  with  the  lower- 
accuracy  fast  exponential  (see  Section  V),  the  argument  must  be  somewhat 
larger.  Since  the  relative  error  approximately  equals  the  argument  of  the 
exponential,  the  criterion  for  using  Eq.  (2.  12)  in  the  PTRANS  subroutine 
is  now  set  at  A  <  2  x  10-2.  With  this  value,  the  greatest  relative  error 
arising  from  the  expansion  and  cancellation  should  be  on  the  order  of 
1  percent. 

2.  1.  3.  Boundary  Conditions 

Integration  of  the  transport  equation  to  obtain  intensities  is  per¬ 
formed  through  the  thickness  of  a  zone,  called  a  "trans"  region.  At 
intersections  of  characteristic  rays  with  the  inner  and  outer  surfaces  of 
each  layer  it  is  necessary  to  supply  the  starting  value  of  the  intensity  I|_i 
required  in  Eq.  (2.  5).  Three  classes  of  boundary  conditions  occur: 


1.  The  trans  region  outside  boundary  coincides  with  outside  zones 
of  the  SPUTTER  calculation  and  a  prescribed  function,  Iq,  is 
applied  at  the  left  boundary  value: 


I(XIA+1’  ~  V 


pL  <  0  or  blackbody  boundary  condition  . 


(2.  13) 
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2.  The  right-hand  boundary  of  the  slab  provide!  for  reflective  and 
transmittal  boundary  conditions  as  well  as  special  routines  to 
establish  prescribed  intensities  for  angles  with  M<0  at  the  boundary 


or 
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I(XIB,  -  p)  =  IOt^,  fi),  i(b1b,  h  )  =  0 
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(2.  14) 


where  Iq  is  the  prescribed  negatively  directed  boundary  intensity 
applied  to  the  outer  ‘bounda.ry  as  a  function  of  angle,  frequency, 
and  time.  Intensities  at  up  to  50  frequencies  and  six  angles  can 
be  accommodated  in  the  table  located1  in  the  array  QINTl(N).  The 
table  entries  are  used  as  Iq  and  are  formed  in  the  subroutine  QUE4 
where  'they  are  stored  in  the  -QINT  1; array.  .  Since  this  subroutine 
is  appropriate  to  the  thermalariteraction -application,  additional 
uses  may  require  subroutines  tailored  to  , the  specific  application. 

AH  other  tians1  boundaries  are  bounded, by  regions  in  which  the 
diffusion  approximation  is  valid  (see  Section  III),  Consequently, 
the  boundary  surface  intensities  ' on  contiguous  trans  regions 
inside  or  outside  of  a  diffusion  region  are  given  by  the  diffusion 
approximation  intensity  derived  in  Section  III: 
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2.  2.  ANGULAR  INTEGRATION 


Integrals  over  the  polar  angle  of  the  intensity  are  required,  as  des¬ 
cribed  in  Section  2.  2. 1,  to  carry  the  calculation  forward  in  time  and  to 
provide  edits  of  informative  derived  quantities.  These  are  formed  by 
numerical  quadrature  using  the  intensities  evaluated  at  a  series  of  discrete 
values  of  polar  angle  by  the  integrations  described  in  Section  2.  2.  Since 
in  the  plane  calculation  ,the  value  of  the  polar  angle  remains  fixed  along  a 
characteristic  ray  and  enters  only  parametrically  in  the  equations,  it  is 
possible  to  exercise  a  choice  of  polar-angle  values  in  order  to  optimize 
the  accuracy  of  the  resulting  integrals. 

The  numerical  quadrature  method  used  for  the  PTRANS  subroutine 
is  the  so-called  double  Gaussian.  ^  In  this  method  the  integrals  of  the 
radiation  quantities  (flux,  energy,  pressure,  etc.  )  are  approximated  by 

f  I  f(fi)  dp  =  V  A  (If)  , 

J  'r '  m  m 

0 
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where  (If)m  is  the  known  value  o£  the  integrand  at  a  chosen  value  pm  of  the 
cosine  of  the  polar  angle,  p.  In  the  method  of  Gaussian,  not  only  are  the 
coefficients  Am  determined  but  the  values  of  pm  are  fixed  to  minimize  the 
difference  between  the  integral  and  the  approximation.  The  result  of  this 
minimization  is  to  relate  the  pm  to  the  zeros  of  the  Legendre  polynomial 
of  order  n  +  1. 

For  those  integrals  having  the  range  -1  <  p  <  1  it  is  frequently 
advantageous  to  treat  the  forward  and  backward  hemispheres  separately 
to  allow  for  the  possibility  of  a  discontinuity  in  I  at  p  =  0..  Such  ’discon¬ 
tinuities  or  very  abrupt  changes  in  the  values  of  the  intensity  between 
forward  and  backward  directions  may  occur  in  systems  which  are  trans¬ 
parent  enough  that  strong  source  regions  are  accessible.  In  these  cases, 
a  better  fit  to  the  integrand  is  obtained  by  the  two  approximating  functions 
which  permit  a  discontinuity  at  p  =  0  than  by  a  single  approximating  function 
which  imposes  a  smooth  behavior  near  p  =  0.  The  method  used  in  PTRANS, 
based  on  separate  integration  regions  for  -1 <  p  <  0  and  0  <  p  <  1,  is  called 
the  double  Gaussian  quadrature  method.  Values  of  Am  and  pm  are  derived 
by  a  simple  transformation  from  those  for  the -single  integration  region. 
Since  the  angle s .for  a  single  integration  region  are  arranged  symmetrically 
about 'the, interval  midpoint,  for  double  region  integration  it  is  possible  to 
identify  pairs  of  angles  ±pm  having  the  same  weight  Am.  In  Tabled,  i,  the 
values(2)  £or  the  0  <.  p  <,  1  interval  are  recorded  for  values  of  n  =  1,  2,t3,  4, 

5.  The  total  number  of  forward  and  backward  angles,  2n  +  2,  for  each  n 
(also  equal  to  the  total  number  of  entries  in  the  table  of  pm  and  Am  for 
each  n)  is  also  listed  in  the  table. 

The  backward  and  the  corresponding  forward  ray  integrations  in  the 
PTRANS  subroutine  are  performed  sequentially.  Since  the  same  absolute 
values  of  pm  are  required  for  these  two  calculations,  many  of  the  quantities 
formed  in  the  backward  integration  pass  can  be  used  for  the  forward  pass 
as  well,  and  hence  these  quantities  are  saved  to  increase  calculation 
efficiency.  Contributions  of  the  pair  of  forward  and  backward  intensities 
to  the  weighted  sums  corresponding  to  the  angular  integrals  are  tallied  at 
the  same  time  that  the  forward  integration  pass  is  being  calculated. 
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SECTION  III 

THE  DIFFUSION  APPROXIMATION 


The  radiation  transport  equation  in  the  limiting  case  of  an  optically 
thick  medium  admits  of  the  diffusion  approximation  in  which  the  expres- 
sion  for  the  radiation  intensity  is  greatly  simplified;  only  the  local  pro¬ 
perties  affect  the  radiation  intensity  at  the  point  in  question.  An  expansion 
of  the  radiation  source  function  Bv  about  the  point  r  permits  the  intensity 
Iy(p)  of  the  radiation  field  in  the  direction  making  an  angle,  whose  cosine 
is  n,  with  the  linear  direction  to  be  formed. 


3.1.  DIFFERENTIAL  FORM  OF.THE  DIFFUSION  ..FLUX 


The  general  solution  of  the  transport  equation  forms  the  starting  . 
point  of  the  derivation;  The  integral  expression  for  the  intensity  appli-. 
cable  to  all  geometries  is 


i(x)  =  f  B(T')e"<T"T,)  dT' 
00 


where  t  =  <vp  ds,  in  which  Kv  is  the  monochromatic. absorption  coef¬ 
ficient  (in  cm^/g)  at  frequency  v.  By  expanding  B(t')  in  series  about  'the 
point  r,  i.  e. , 

2 

,i  „/  »  .  9B  ,  ,  ,  .  I  B  B,  .  .2  . 

B(t')  =  B(r) +  ^- ^t'  -  t) +  - — j-(t'-t)  +  •••  , 


the  intensity  becomes 


9t 


T-B  ®£  +  iI-2 

1-8  28r2 


or 


=  B .  jl*s +il±  (IlmV 

Kp  9x  Kp  8x  \Kp  9x  J 


for  plane  slab  geometry. 

The  diffusion  approximation  results  from  retention  of  only  the  first 
two  terms,  so  that  the  diffusion  intensity  is 


(3.1) 


I  =  B 


8B 
<p  8x 


and  the  monochromatic  diffusion  flux  ^  and  radiation  energy  E  are 

r  xv 


4> 


r 


4trc  1  8B 
3  Kp  8x 


(3.2) 


3.  2.  CRITERIA  FOR  THE  SELECTION  OF  DIFFUSION  REGIONS 


The  criteria  for  the  validity  of  the  diffusion  approximation  can  be 
obtained  by  examination  of  the  above  derivation- -namely,  that  the  expan¬ 
sion  of  the  source  function  be  justified- and  that  the  expansion  converge 
rapidly  so  that  the  neglect  of  all. but' the  Leading  terms  is  valid.  If  the 
source  function  is  linear  in  t*  at  the  point  in  question  and  is  also  linear 
for  a  distance  of  the  order  of  one  mean  free  path  on  either  side  of  the  point, 
the  criteria  are  satisfied.  These,  criteria  are  difficult  to  quantify  since 
they  refer  to  a  finite  region  containing  the  point  in  question.  If  all  of  the 
terms  (or  a  large  number  of  them)  were  checked  for  rapid  convergence, 
this  would  imply  (making  a  smoothness  assumption)  that  the  diffusion 
criterion  is  met.  It  is  not  possible  with  finite  differences,  however,  to 
form  the  higher-order  local  derivatives  approxixnationa. 

In  the  SPUTTER  subroutine  PRADTN,  criteria  designed  to  give  an 
indication  of  both  the  local  and  nonlocal  behavior  h&ve  been  employed. 

First,  at  the  zone  interface  at  which  the  intensity  and  flux  are  to  be 
evaluated,  the  inequality 


is  required.  In  this  expression  h  =Sicp  dx  is  the  optical  depth  normal  to 
the  slab;  the  derivative  is  approximated  by  the  centered  first  difference  of 
B  between  adjacent  zones.  The  resulting  expression,  of  course,  contains 
some  nonlocal  aspects  resulting  from  the  finite  difference  approximation, 
which  ensures  that  when  neighboring  zones  are  optically  thick,  no  nonlocal 
source  perturbation  is  close  enough  to  invalidate  the  diffusion  approximation. 
However,  to  provide  for  the  cases  when  a  source  perturbation  is  located 
a  fraction  of  an  optical  depth  from  an  interface  meeting  the  condition  of 
Eq.  (3.  3),  the  diffusion  region  is  constricted.  Starting  froftfthe  closest 
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interfaces  outside  the  diffusion  region  (where  Eq.  (3.  3)  is  not  satisfied), 
all  of  those  interfaces  lying  within  a  prescribed  number  of  mean  free  paths 
are  removed  from  the  diffusion  region. 


The  criteria  used  in  SPUTTER  are  controlled  by  input  numbers. 
The  criterion  of  Eq.  (3.  3)  uses  the  input  number  HCB: 


jTG  |  <  HCB  x  Y2  , 


(3.  4) 


where  TG  is  the  difference  approximation  to  the  gradient  and  Y2  is  the 
source  function  evaluated  at  the  interface  by  interpolation.  The  second 
criterion  uses  the  input  number  HVB  (in  mean  free  paths).  If 


|Q3^  )  -  Q3(J) j  >  HVB  , 


(3.5) 


then  the  interface  with  index  J  which  satisfies  Eq.  (3.  4)  is  removed  from 
the' diffusion  region.  In  Eq.  (3.  5);  Q3  is  the  normal  optical'  depth  and  I  is 
the  index  of  the  nondiffusion  interface  adjoining  the  diffusion  region. 

Although  the  diffusion  calculation  is  considerably  faster  than,  the 
transport,  the  establishment  of  two  transport  regions  separated  by  the 
single  zone  requires  still  more  calculation  to  set  up  characteristic  rays 
and  perform  bookkeeping  operations.  To  avoid  the  duplicate  setup  calcu¬ 
lations  required  for  an  additional  transport  region,  a  test  is  made  to 
eliminate  a  diffusion  region  consisting  of  a  single  zone. 


3.  3.  DIFFERENCE  FORM  OF  THE  DIFFUSION  FLUX 

The  diffusion  intensity  derived  above  is 

\  •  *•  \ 

.  u  3B 

I  =  B  -  -£— -r —  .  <■ 

icp  9x 

In  the  group  frequency  approximation  of  SPUTTER,  the  intensity  integrated 
over  a  frequency  interval  (vj,  v^+j)  is  required: 

Jv.  Jv.  -'v.  99  v 

3  3  3 


In  terms  of  the  partial  Rosseland  mean  absorption  coefficient 
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the  frequency  group  intensity  becomes 


v 


It  is  desired  to  evaluate  this  quantity  at  each  zonal  interfacein  the, 
mesh.  Since  Lhe  known  quantities  are  the  sore  temperatures  and  densities, 
the  absorption  coefficients  Kj  and  the  integrated  source  functions  X6  =S  Bdv 
are  first  evaluated,  not  at  the  interfaces  but  at  positions  representative  of 
each  zone. 

The  question  remains  as  to  how  best  to  approximate  the  derivatives 
and  interpolate  for  the  coefficients  in  Eq.  (3.  6)  at  the  interfaces  from  the 
quantities  available  at  zone  positions.  The  answer  depends  on  the  tem¬ 
perature  and  density  profile  across  the  interface  from  which  these  terms 
could  be  calculated  directly.  Since  the  profile  is  not  known,  we  must 
select  a  reasonable  approximation  which  will  permit  the  calculation  to  be 
carried  out.  In  fact,  the  appropriate  profile  depends  on  the  events  which 
have  taken  place  in  the  calculation  and  on  the  energy  transport  mechanisms 
of  greatest  importance  in  it.  As  extreme  examples,  a  problem  dominated 
by  hydrodynamics  might  have  quantities  determined  by  passage  of  a  strong 
shock  and  subsequent  linearization  in  mass  coordinates  of  the  pressure 
behind  the  shock,  whereas  a  radiation- dominated  diffusion  problem  is 
characterized  by  linearity  of  the  radiation  potential,  which,  in  turn, 
depends  on  the  Rosseland  opacity.  Of  course,  such  detailed  information 
about  the  progress  of  a  problem  is  generally  unavailable,  so,  at  best,  an 
approximation  based  on  over-all  accuracy  is  needed. 

Since  the  terms  under  consideration  are  the  radiation  diffusion 
equations,  the  interpolation  is  performed  in  a  way  to  give  greatest  accuracy 
when  the  diffusion  terms  are  most  important- -namely,  when  the  profile  is 
being  determined  entirely  by  radiation  diffusion.  It  is  also  desirable  to 
reduce  the  number  of  coefficients  requiring  interpolation.  This  can  be 
done  by  noting  the  identity 


and  by  forming  the  variable  t  =  Jp/Cj  dx.  In  terms  of  these  quantities,  the 
intensity  can  be  written  as 
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SECTION  IV 

FREQUENCY  INTEGRATION 


Equations  derived  in  Sections  II  and  III  which  are  applicable  to  a 
particular  frequency  of  the  radiation  field  are  of  limited  usefulness  in  the 
SPUTTER  calculations.  Although  in  principle  a  calculation  at  a  particular 
frequency  might  be  valuable  for  comparison  with  high-resolution  spectro¬ 
scopy,  in  practice  no  such  data  have  been  available.  Of  much  more  use 
are  intensity  s  averaged  over  a  wide  frequency  band.  These  quantities  can 
be  compared  with  data  from  wide-band  measurements  and,  most  important 
of  all,  can  be  summed  for  use  in  the  energy  integration  in  the  SPUTTER 
code.  The  quantities  to  be  summed  are  the  frequency-integrated  radial 
flux  component,  the  radiation  energy  density,  and  the  radiation  pressure. 
For  performing  interaction  calculations,  it  is  also  valuable  to  form  other 
components  of  the  radiation  flux. 

Basically,  the  quantity  which  is  required  for  each  of  the  above 
applications  is  the  frequency-group  intensity  I. ., 


(4.1) 


Then,  for  example,  this  quantity  can  be  integrated  over  angles  to  form 
4>. .,  the  contribution  to  the  flux  at  position  i  of  frequency  group  j: 


*  ■  /: 


1..)  M  dM 


and  thus  the  total  radiant  flux  at  position  i  is 


4>.  =  5*  <J>. .  . 

1  .  1J 

J 

Equation  (2.  8)  gives  the  expression  for  the  frequency-dependent  intensity 
to  be  used  in  Eq.  (4.  1).  ^he  frequency  integ.'  ♦'ion  of  Eq.  (2.  8)  has  been 
reported  recently,  (4)  but  the  current  SPUTTExv  code  does  not  include  the 
transmission  functions.  The  first  two  terms  of  Eq.  (2.  8)  which  form  the 
diffusion  limit  can  be  integrated,  as  in  Section  3.  3,  to  give 
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{diffusion  limit)  , 


(4.2) 


i.: 

ij 


3.. 

ij 


o-  9x 
R. 
j 


in  which  the  first  term 


B.. 


B^(v)  dv 


is  the  frequency- group  Planck  function  and  the  second  term  contains  the 
frequency-group  Rosseland  mean  absorption  coefficient  <tr.  -  9  <y  In  this 
form,  Eq.  (4.  2)  correctly  gives  the  frequency-group  intensity  for  the 
optically  thick  limiting  case.  The  remaining  Bj  and  8B/8T).  terms  of 
Eq.  (2.  8)  are  formed  in  the  same  way.  Thus, 


In  Eq.  (4.  3),  mean  values  of  the  exponentials  have  been  extracted  from  the 
frequency  integrals  and  the  outstanding  problem  is  to  specify  their  values. 
Two  options  are  available;  they  differ  in  the  absorption  coefficient  used  to 
calculate  the  optical  depth.  The  first  is 


-A 

e 


and  the  second  is 


(4.  4) 


where 


and 


f  lj+1  cr  B  dv 
J  V  V 

v  . 

_J _ 

B.. 


6  = 


x. 

l 


-  X 
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For  small  optical  depth,  the  correct  result  makes  use  of  the  Planck  mean 
absorption  coefficient.  From  Eq.  (2.  12)  the  frequency  integration  then 
gives 


I.. 

lJ 


=  I. 


i-1.  j 


4b..  +4b.  .  .  +  k  ,  , 

_4  1J  4  1-1,  j  2  i-£,j 


(4.  5) 


The  above  prescriptions  for  frequency-group  means  are  far  from 
satisfying  and  call  for  further  work.  Considerable  economies  can  be  made 
through  reductions  in  the  number  of  frequency  groups  if  a  more  accurate 
means  of  averaging  within  groups  can  be  found.  Presently  used  choices 
of  frequency  groups  appear  to  give  a  reasonably  accurate  result,  however, 
as  indicated  by  comparisons  between  calculations  with  the  nominal  number 
of  frequency  groups  and  calculations  with  a  very  large  number  of  frequency 
groups.  (It  is  expected  that  a  unique  correct  result  will  be  obtained  as  the 
number  of  frequency  groups  is  increased,  irrespective  of  the  choice  of  the 
weighting  function  in  the  frequency-group-average  absorption  coefficient.  ) 
Consequently,  a  very  few  frequency  groups  should  be  adequate  if  a  suitable 
averaging  procedure  were  developed. 


Even  with  a  crude  averaging  scheme,  considerable  improvement  in 
accuracy  results  from  choice  of  frequency- group  boundaries  so  as  to 
reduce  the  variation  of  the  absorption  coefficient  within  the  group. 


Work  on  the  absorption  coefficient  for  air  indicates  that  approxi¬ 
mately  20  groups,  carefully  selected  as  to  their  locations,  afford  quite 
adequate  resolution.  Enough  information  is  known  about  air  to  make  this 
selection  appear  quite  reasonable.  Air  absorption  coefficient  tapes 
(DIANE)*  have  been  prepared  for  18,  20,  and  90  groups.  The  90-group 
tape  is  used  to  check  on  the  frequency  integrations  at  selected  times.  The 
proper  averages  to  use  are  difficult  to  decide  on  at  this  time.  There  are 
provisions  for  reading  into  storage  from  the  DIANE  tapes  both  the  Rosseland 
and  Planck  averages,  which  are  used  at  present  in  the  thick  or  thin  limits, 
respectively. 


See  Section  VI  of  Volume  V. 
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SECTION  V 

SUBROUTINE  ORGANIZATION  AND  ECONOMICS 


The  present  plane  transport  subroutines  were  written  with  the  idea 
of  removing  unnecessary  calculations  from  inside  the  frequency  loop  and 
characteristic  ray  integrations.  These  improvements  required  am  increase 
in  storage  for  the  subroutines  to  attain  a  decrease  in  calculational  time. 

The  reorganized  subroutines  will  be  discussed  in  two  sections,  corres¬ 
ponding  to  the  two  major  subroutines:  (1)  the  radiation  subroutine  (PRADTN) 
in  which  most  of  the  preliminary  setup  and  the  diffusion  calculation  is  com¬ 
pleted  and  (2)  the  transport  subroutine  (PTRANS)  in  which  the  intensity 
calculation  and  angular  integrations  are  performed.  The  subroutines 
which  execute  the  opacity  interpolations  (KAPPA),  Planck  function 
(PLNKUT),  and  fast  exponential  (FREXP)  will  be  discussed  in  Section  5.4. 
The  input  numbers  and  the  output  edits  will  be  presented  in  Sections  5.  5 
and  5.  6. 


5.  1.  THE  PRADTN  SUBROUTINE 

In  PRADTN,  the  high-frequency  groups  are  merged,  a  source  region 
is  established,  boundary  sources  and  derivatives  are  calculated,  regions 
for  transport  and  diffusion  are  formed,  diffusion  fluxes  are  calculated, 
frequency  integration  is  performed,  and  the  radiation  time-step  control  is 
evaluated.  Each  of  these  activities  in  PRADTN  will  be  discussed  in  sub¬ 
sequent  paragraphs. 

5.  1.  1.  Merge  Frequency  Groups 

Frequency  groups  that  are  too  far  out  on  the  Planck  tail  for  a 
"maximum"  temperature  in  the  mesh  are  merged.  The  criterion  used  is 
as  follows:  If  the  lower  frequency  boundary  hvj  of  the  group  in  question 
(hvj,  hv£)  is  greater  than  ten  times  the  maximum  temperature  (THMAX) 
in  the  mesh,  this  group  will  be  merged  with  the  next  lower  group.  Merging 
will  continue  until  over  half  the  groups  have  been  merged;  at  this  point, 
either  the  calculation  is  terminated  or  a  second  DIANE  tape  is  called.  On 
merging,  Rosseland  and  Planck  averages  are  formed  by  using  the  following 
equation  for  dB/d0^  and  the  appropriate  sums: 
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4  dB  4  dB 

I  b;e  *  I  — r.  I  v  S’  I  — r~  •  t5-1) 

J  d0  J  dQ  Kr 

The  Planck  weighting  functions  (bj)  are  obtained  from  PLNKUT,  as  des¬ 
cribed  later.  On  completing  the  merging,  the  merged  opacities  are 
formed: 


>cR  =  £  dBy/d e4/ £  dBv/(de4 


(CAPAR) 


(5.2) 

(CAPAC)  . 


The  frequency- dependent  sources  must  be  established  at  the  inter¬ 
faces  from  the  zonal  quantities  bj0j*+i  (X6(i))  and  t^I  (H3(i)).  The  dif¬ 
ference  equations  used  were  given  in2  Section  2.  1.  Before  the  calculation 
of  the  Planck  function  (bj)  is  made,  i.  e. ,  before  calling  PLNKUT,  a  test 
is  made  to  see  if  uj  (i.  e.  ,  the  reduced  frequency  hvj/0)  >  19;  if  so,  bj  =  0 
(i.  e. ,  the  source  X6(i)  =  0.  0).  If  uj  <  19.  0  and  u;>  <  0„  01,  then  b j  =  0  also, 
assuming  that  for  0*  <  105,  the  small  bj  (bj  ~  10"5)  will  produce  a  negli¬ 
gibly  small  source  contribution.  An  index  (ICX)  is  set  equal  to  the  last 
zone  that  contains  a  source.  This  source  index  is  used  to  limit  the  trans¬ 
port  calculation  to  the  region  containing  sources.  While  setting  up  the 
sources  and  derivatives,  tests  are  made  on  their  discontinuous  nature  to 
use  either  a  linear  or  constant  form  in  the  intensity  integrations.  7ihe 
initial  check  is  on  the  minimum  optical  depth  of  adjacent  zones  to  ensure 
that  both  are  transparent  (less  than  0.  3).  If  this  condition  holds  and  if 
both  the  sources  and  optical  depths  are  changing  rapidly  in  x  (change 
greater  than  a  factor  of  two),  the  derivative  at  that  interface  (TG(i))  is  set 
equal  to  zero.  The  zero  source  derivative  is  used  in  PTRANS,  as  a  test, 
to  set  up  the  constant  source  terms.  For  the  intensity  integration,  special 
boundary  sources  and  derivatives  are  also  established  at  the  edge  of  the 
source  region  (I  =  ICX)  and  at  the  outside  of  the  mesh  (I  =  IM)  (see 
Section  2.  1 .  3). 


5.1.3.  Determine  Diffusion  Region 

The  principal  criterion  for  defining  a  diffusion  region  is  that  the  first 
derivative  of  the  source  function  (TG)  be  small  compared  to  the  source 
(Y2)  {see  Section  3.  2).  When  the  zone  is  found  to  be  diffusion,  the  boundary 
is  tagged  by  setting  (X3  =  -1. ).  Before  incorporating  this  interface  into  a 
diffusion  region,  the  possible  influence  from  sources  on  either  side  is  con¬ 
sidered  and  a  further  test  is  made.  From  the  last  diffusion  boundary,  a 
test  is  made  for  an  optical  depth  in  succeeding  zones  to  the  left.  If  more 
than  HVB  optical  depths  appear  in  the  next  zone,  then  this  zone  is  calculated 
by  transport  and  removed  from  the  diffusion  region  (set  X4(i)  =  -1.  0). 

HVB  is  an  input  number,  which  is  usually  around  5.  When  x  =  0  is  reached 
after  testing  each  zone,  zones  out  to  the  right  of  the  present  transport 
region  are  tested  in  the  same  manner.  The  above  test  buffers  the  trans¬ 
port  region  with  an  (HVB)  mean-free- path- thick  diffusion  boundary.  If  the 
zone  boundary  stays  diffusion,  i.  e. ,  X3(i)  =  -1.0  and  X4(i)  =  0.  0,  a  dif¬ 
fusion  flux  is  calculated  from  the  source  gradients,  as  described  in 
Section  3.  1.  The  regions  where  X3(i)  =  0.  or  X3(i)  =  -1.  and  X4(i)  =  -1. 
have  been  established  as  transport  regions  because  they  did  not  meet  the 
diffusion  criteria  or  they  reverted  to  transport  regions  by  the  optical- 
depth  test  described  above.  This  transport  region  is  then  identified  by 
setting  the  left  boundary  to  LAX  and  the  right  boundary  to  IBX.  More  than 
one  trans  region  may  be  set  up  in  PRADTN,  and  if  so,  a  PTRANS  calcu¬ 
lation  will  be  made  for  each  region.  No  one-zone  diffusion  region  is 
allowed  and  the  region  outside  the  sources  (I  >  ICX)  is  always  considered  a 
transport  region. 

5.  1.4.  Time -step  Control  and  Monofrequency  Calculation 

These  two  aspects  of  the  new  code  are  related  since  the  "grey" 
absorption  coefficients  from  the  DIANE  tape  are  used  to  estimate  a  radi¬ 
ation  time  step  as  well  as  to  form  the  monofrequency  time -dependent 
calculation.  In  the  multifrequency  calculation,  after  all  groups  have  been 
processed,  an  additional  call  for  KAPPA  is  made  to  read  in  the  grey 
absorption  coefficients.  These  averages  were  obtained  by  integrating  the 
frequency- dependent  absorption  coefficients  for  both  Planck  (/Cp)  and 
Rosseland  (k%)  in  the  DIANE  code.  The  actual  time  step  for  radiation 
transfer  is  then  obtained  from  the  formula 

At  =  (0.  5  +  1.  5  H3(i)2)/(acK_03)x  CV(i)  ,  (5.3) 

K 

12 

where  CV(i)  is  the  specific  heat  and  ac  =  4.  12  x  10  .  The  mass  point  in 

question  is  also  checked  to  ensure  that  it  will  not  gain  or  lose  more  than 
half  its  original  energy: 

At-  =  0.  5  x  CV(i)  x  0(i)  x  G(i)/|ER(i)|  ,  (5.4) 

K 
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where  £R(i)  is  the  divergence  of  the  flux  and  G(i)  is  the  mass  in  the  zone. 
The  minimum  of  these  values  is  compared  to  the  hydro  time  step  (Courant) 
and  if  smaller, 

NRAD  =  FIX(DTH2/DTRMIN)  and  DTR  =  DTH2/NRAD  (5. 5) 
is  set  to  cycle  NRAD  times  through  the  radiation  routine. 

The  monofrequency  calculation  also  uses  the  grey  absorption  coef¬ 
ficients  from  the  DLANE  tape.  If  KMAX  =  0.  0  and  S15  =  1.0,  the  frequency- 
averaged  opacities  are  bypassed  on  the  tape  and  only  the  grey  absorption 
coefficients  are'readinto  storage.  For  succeeding  cycles,  S15  is  set 
equal  to  zero  and  the  interpolations  for  and  Kp  are  performed  in 
KAPPA  using  the  stored  opacities  originally  read  into  KAPPA's  common 
storage.  When  the  problem  is  restarted  it  is  therefore  necessary  to  reload 
SI  5  equal  to  one.  If  the  DLAN£  tape  is  not  designated  (the  tape  unit  assigned 
must  be  stored  in  AMASNO(J+l  7),  where  J  is  the  material  number),  then 
the  KAP  routine  is  called  (KAP8  for  air)  and  used  for  the  monofrequency 
calculation. 


5.2.  THE  PTRANS  SUBROUTINE 

The  subroutine  PTRANS  is  called  by  PRADTN  to  carry  out  the 
intensity  integration  between  LAX  and  IBX,  saving  various  quantities  on 
the  inward  pass  that  will  be  used  on  the  outward  pass  as  well  as  the 
angular  integration  of  the  flux  between  rays  ( J'.}  If*  dji).  After  the  inten¬ 
sity  transport  along  a  typical  ray  in  the  outward  direction  (iA  -*■  iB)  is  done, 
the  flux  is  calculated  while  the  inward  pass  of  the  intensity  calculation  is 
being  completed.  The  angular  integration  is  based  on  a  linear  inter¬ 
polation  of  the  intensities  between  rays.  The  logic  in  PTRANS  is  des¬ 
cribed  in  detail  in  the  following  sections. 

5.  2.  1.  Selection  of  Angles 

At  present,  only  five  sets  of  Gaussian  angles  and  weights  are  stored 
in  the  subroutine.  These  can  be  selected  by  setting  an  input  number 
(LMDA(37),  the  number  of  angles  with  fi> 0)  to  the  desired  n+1.  The 
selection  from  storage  is  made  from  the  following  indices 

NY  =  LMDA(37)  -  1 
NMU  =  (NY-1)  x  (NY+2)  +  1 
NGS  =  NMU+NY  +  1  , 

NMU  selects  the  cosine  of  the  angle  (Mm);  NGS  selects  the  relation 
(^m  Am)»  the  cosine  of  the  angle  times  the  Gaussian  weights  for  the  flux 
formulation. 
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5.  2.  2.  Intensity  Integration  along  Characteristic  Rays 


The  integration  using  Eq.  (2. 11)  starts  at  the  left  boundary  with  the 
appropriate  boundary  condition  and  proceeds  outward,  storing  the  expo¬ 
nentials  e"^T  in  (H4(i)),  the  derivatives  \i  8B/9h)i+IinX8(i),  and  the  cal¬ 
culated  intensities  in  sum  X3(i).  The  more  general  boundary  conditions 
are  established  (see  Section  2.  1)  and  the  stored  quantities  are  now  used 
except  for  the  change  of  sign  of  n  8B/8h)i+i  in  Eq.  (2.  1)  to  calculate  the 
intensities  I(F2),  on  the  outward  pass. 

The  regions  where  constant  sources,  and  therefore  zero  boundary 
derivatives,  should  be  used  in  the  intensity  integrations  were  established 
in  PRADTN  by  setting  TG(i)  equal  to  zero.  In  the  integration  along  a 
particular  ray,  a  test  is  made  on  TG{i)  at  each  interface;  if  zero,  the 
source  terms  Y2(i  -  1)  and  Y2(i)  are  set  equal  to  X6(i  -  7)  respectively 
(see  Fig.  2.  1). 

As  discussed  in  Section  2.  1.  1,  the  accuracy  of  the  exponential  term 
and  the  effect  of  truncating  errors  mean  that  the  general  formula  will  not 
reduce  in  the  limit  of  small  optical  depths  to  the  transparent  case.  To 
correct  this  situation,  a  test  is  made  on  tj  (the  half  optical  depth  t  is 
stored  in  H2(i)),  and  if  r  <  1 0”2  a  switch  is  made  to  the  limiting  form  of 
the  transport  equation  (Eq.  (2.  12)  developed  in  Section  II). 

5.  2.  3.  Angular  Integration 

The  only  integral  over  angle  formed  in  PTRANS,  at  present,  is  the 
flux;  (if  In  dp)  the  formula  for  energy  (vfldn)  is  included  for  possible  use 
later.  These  integrals  are  formed  on  the  outward  pass  from  the  intensity 
(sum  X3(i))  stored  on  the  inward  pass  and  the  intensity  being  calculated 
(F2).  The  difference  forms  of  the  equations  are 

X2(i)  =  V  (F2  -  sum  X3(i))  x  u  A  , 

*—  mm 

ER(i)  =  I  (F2  +  sum  X3(i))  x  A 

m 

5.  3.  DIFFERENCES  WITH  INTEGRAL  FORMULATION 


The  principal  difference  in  the  subroutines  is  in  replacing  the  inte¬ 
gration  of  angle  done  explicitly  in  the  integral  formulation  by  a  sampling 
scheme  of  a  double  Gaussian  nature.  It  is  expected  that  accuracy  can  be 
achieved  with  a  minimum  number  of  rays  (presumably  less  than  n  =  6,  see 
Section  2.  2).  This  result  is  in  logical  argeement  with  the  use  of  the  S4 
approximation  in  the  neutron-transport  work.  The  advantage,  therefore, 
will  appear  in  problems  with  many  zones,  since  the  integral  method  will 


increase  as  the  number  of  zones  squared  whereas  the  present  method  will 
only  increase  linearly  with  zones.  Furthermore,  the  present  method 
makes  it  possible  to  have  special  boundary  conditions  depending  on  angle 
(see  Section  2. 1.  3). 


5.  4.  AUXILIARY  SUBROUTINES 


In  addition  to  the  two  new  basic  subroutines  PRADTN  and  PTRANS, 
some  changes  have  been  made  in  the  auxiliary  subroutines  EXP,  PLNKUT, 
and  KAPPA.  These  changes  include  (1)  a  fast  exponential  (FREXP),  (2)  a 
two-argument  Planck  function,  and  (3)  the  use  of  the  average  opacities 
from  KAPPA  (6  and  p  interpolation)  for  the  monofrequency  calculation  as 
well  as  for  the  Planck  opacities. 

The  new  fast  exponential  routine  FREXP  uses  table  lookup  and  inter¬ 
polation  rather  than  the  normal  expansion  methods.  The  routine  is  written 
in  machine  language  but  uses  the  library  routine  EXP(X)  for  positive  X  or 
X  >  -10.  An  over-all  gain  in  speed  of  a  few  percent  was  achieved  in  one 
comparison  SPUTTER  calculation. 

The  PLNKUT  routine,  with  its  associated  tables  PLNKTT,  has  been 
corrected  and  made  more  efficient  by  using  a  two-argument  call  which 
now  calculates  from  either  the  analytic  form  or  from  the  tables  the  dif¬ 
ference  in 


b<ur  f 2  *£■--- C7Er-7*’-bj  •  (5-6) 

Jv j  c  e  -  1 

The  accuracy  is  improved  since  now  not  only  differences  of  nearly  equal 
numbers  are  subtracted. 


The  subroutine  KAPPA,  which  calls  in  the  group-averaged  absorption 
coefficients  from  the  DIANE  tape  and  performs  a  bilinear  log  interpolation 
in  temperature  and  density,  has  been  modified  to  obtain  the  grey  absorption 
coefficients  as  well  as  the  Planck  averages.  At  present,  the  format  of 
the  DIANE  (absorption  coefficient)  tape  includes  a  BCD  record  for  tape 
identification,  the  Rosseland  and  Planck  averages  for  a  selected  set  of 
temperatures  and  densities  from  0.  25  ev  to  50  ev  and  from  10  normal  to 
10-6  normal,  and  the  actual  integration,  J*  Kv  dv,  for  the  grey  case.  The 
grey  or  frequency-integrated  averages  are  also  used  for  an  estimate  of 
the  time  steps  in  PRADTN.  KAPPA  reads  in  first  the  tape  name,  the 
number  of  frequency  groups,  and  the  size  of  the  records.  If  the  sentinel 
for  multifrequency  is  set  to  KMAX  =  1,  then  the  first  frequency  group, 
hvi?  and  its  absorption  coefficients  are  read  into  storage,  The  inter¬ 
polations  in  log  ©i  and  log  are  performed  and  a  return  to  PRADTN  is 
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made.  If  KMAX  =  0,  then  KAPPA  skips  over  the  frequency- dependent 
absorption  coefficients  and  reads  into  storage  the  grey  averages.  A 
signal,  SI 5  =0,  is  subsequently  set,  and  for  further  cycles  the  inter¬ 
polations  are  made  on  the  stored  quantities;  the  tape  is  not  called  again. 

5.  5-  INPUT  NUMBERS 

The  input  quantities  used  in  the  radiation-transport  subroutines  and 
their  functions  are  listed  in  Table  5.  1.  The  entries  in  it  are  as  follows: 
column  1  is  the  storage  location  number  used  for  entering  the  quantity  into 
storage  with  the  CARDS  subroutine,  column  2  lists  the  FORTRAN  name  of 
the  stored  quantity,  column  3  gives  the  range  of  admissible  values  of  the 
input  quantity,  column  4  describes  its  function  and  identifies  special  values 
it  may  assume,  and  column  6  records  a  set  of  values  of  the  quantities  which 
might  be  typical  of  those  for  a  normal  problem.  Included  is  a  set  of  values 
for  the  input  quantities  selected  for  solving  typical  problems. 

5.  6.  EDITS 

The  editing  of  such  frequency- dependent  quantities  as  H3,  the  optical 
depth  (Rosseland),  X6,  the  source  (bju  )»  X2,  the  flux  (in  ergs/4/3  it  sec) 
X2/DKNU,  the  flux  divided  by  the  frequency  group,  THETA,  the  temper¬ 
ature  (ev),  and  El,  the  energy  (ERGS/G)  versus  radius  is  accomplished  by 
setting  S12  to  the  desired  number  of  cycles  between  prints.  These  multi¬ 
frequency  edits  have  been  used  to  evaluate  the  criteria  for  the  subroutines 
as  well  as  for  diagnostics  during  the  calculations. 

A  list  of  sample  editing  for  a  particular  frequency  group  is  given  on 
page  27.  The  HNU  is  in  electron  volts.  The  quantities  found  useful  to  dis¬ 
play  for  each  frequency  group  and  for  a  characteristic  ray  are  listed  on 
page  28.  The  format  statements,  in  the  listings  appended,  have  been 
revised  for  the  debug  print  from  those  used  on  page  28. 
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Table  5.  1 

PLANE  RADIATION  INPUT  QUANTITIES 


Card 

No. 

Quantity 

Range 

of 

Values 

1 

Description 

Typical 

Value 

37 

LMDA(37) 

2.  3,  4.  5.  6 

Number  of  angles  with  ft  >  0. 

2 

44 

KMAX 

#  0,  performs  multigroup  frequency  approximation; 

=  0,  performs  single-group  frequency  approximation. 

0 

81 

HVB 

>  0 

Number  of  optical  mean  free  paths  by  which  a  transport 
region  is  extended  at  the  expense  of  each  adjacent  diffusion 
region.  (See  Section  3.  2.  ) 

5 

83 

■ 

HCB 

>  0 

Criterion  to  define  a  diffusion  region  in  terms  of  relative 
gradient  of  the  source  function.  Diffusion  regions  are 
eliminated  if  0.  (See  Section  3. 2. ) 

0. 1 

87 

CB 

>  0 

Criterion  to  combine  frequency  groups.  If  the  lower  fre¬ 
quency  of  the  group  is  more  than  CB  times  the  temperature 
of  the  hottest  sons,  that  group  is  combined  with  the  adjacent 
greup  of  lower  frequency.  A  half-integer  value  presents 
termination  of  the  problem  when  half  or  more  of  the  groups 
have  been  combined.  (See  Section  5. 1. 1. ) 

10.5 

88 

GA 

£  0 

One  of  two  criteria  for  choice  of  linear  or  stepwise  constant 
source  within  a  zone.  (See  Section  5. 1.2.) 

0.  3.33 

90 

GL 

Neg. ,  0, 

i.  p°»- 

integer 

Indicator  for  radiation  boundary  condition  at  IB. 

GL  =  negative,  total  reflection; 

GL  =  0,  intensity  for  p  <  0  is  zero; 

GL  =  blackbody  intensity  based  on  temperature  located 

in  THETA(IB)  for  m  <  0; 

GL  =  positive  integer,  intensity  for  p  <  0  obtained  from 

source  routine.  GL  must  equal  number  of  frequency 
groups. 

(See  Section  2.  1.  3. ) 

0 

121 

AC 

>  0 

One  of  two  criteria  for  choice  of  linear  or  stepwise  con¬ 
stant  source.  Minimum  value  for  using  a  linear  source. 

(See  Section  5.  1.  2. ) 

0.  3 

127 

AC03T4 

— 

Transport  debug  edit  criterion.  Edit  occurs  if  i  0  and 
<  cycle  number.  | 

9 

147 

S12 

---- 

Number  of  cycles  between  multifrequency  edits. 

10 

150 

S15 

0.  t 

Trigger  controlling  call  of  DIANE  tape.  Must  have  value 
!  d  0  on  starts  or  restarts. 

1 

8466 

T£JLM(25) 

Constant  multiplying  the  radiation  time  step.  Can  be  used 
to  modify  the  stability  criterion. 

1 

8478 

TELM(37) 

>0 

Maximum  permissible  fractional  energy  in  any  zone  due  to 
radiation.  Time  step  may  be  reduced  to  meet  this 
requirement. 

0.  05 

8858 

SOLID(IO) 

! 

Thick-thin  criterion,  If  0,  Planck  mean  is  used  to  form 

H2;  otherwise,  Rosseland  mean.  (See  Section  4.  ) 

1 
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6622066*01  9.23704*5-02 


Table  5.  2 

SAMPLE  MULTIFREOUENCY  EDIT  FOR  PTRANS 
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SECTION  VI 


TIMING  STUDIES  AND  ACCURACY  IN 
ANGULAR  INTEGRATION 


6.  1.  TIMING  CALLS 

The  comparison  in  calcul  don  time  wai  obtained  by  using  timing 
calls  at  selected  locations  in  the  logic  of  the  code.  To  use  the  timing  calls, 
it  was  necessary  to  establish  a  fiducial  time  from  the  system  clock  and 
then  print  the  location  of  the  time  call,  the  time,  and  the  difference  in 
time  between  calls  for  each  call.  The  subroutine  that  carries  out  these 
steps  is  CLOCK. 

In  the  calculations  described  above,  the  subroutine  CLOCK  was 
called  at  the  following  locations  in  PRADTN  and  PTRANS: 


PRADTN 


13. 

105 

13. 

140 

13. 

701 

13. 

151 

13. 

180 

13. 

292 

13. 

286 

-13. 

239 

Before  frequency  loop 

After  call  KAPPA  on  merge 

After  call  KAPPA  on  main  frequency  loop 

After  calculating  general  sources 

Before  calling  PTRANS 

After  EDIT  (normal)  end  of  frequency  loop 

After  last  frequency  start  time  step 

End  of  cycle  (return  to  main  program) 


PTRANS 


14.  708  -  Before  debug  print 

The  following  calculation  was  timed  in  units  of  1/60  sec  for  the 
above  breakdown  in  computing  time. 


The  calculation  described  here  did  not  use  the  TG  criteria  (see 
Section  5.  1. 2)  nor  the  special  boundary  conditions  (see  Section  5.  2.  4). 
Three  ray  passes  were  completed  for  each  frequency  group  for  a  total  of 
six  angles,  forward  and  back,  and  with  32  active  zones.  The  total  time 
for  21  frequency  groups  was  ~14.  9  sec.  The  breakdown  in  time  for  a 
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single  frequency  group  (in  units  of  1/60  sec)  are  the  following: 


Call  KAPPA  for  absorption  coefficient 
After  sources 

Characteristic  ray  passes  (3) 

After  call  PTRANS  with  EDIT 
Average  time  required 


uX 

1. 

3. 

23. 

'41 /hv 


\ 


The  start  and  merge  of  KAPPA  is 
Total  time  with  EDIT 
Total  time  without  EDIT 


~32 

~869  (1/60  sec) 
~51 5  (1/60  sec) 


X 


X 


\v 


\ 


6.2.  ACCURACY  IN  ANGULAR  INTEGRATIONS 

Comparisons  have  been  made  between  calculations  for  two  sets  of 
angles  for  a  radiation  shock  problem.  The  problem  consists  cl  a  hot 
(5  ev)  shock  moving  into  cold  low-density  air.  The  effect  on  the  flux 
versus  linear  zones  for  a  set  of  four  and  eight  angles  is  given  in  Fig.  6.  1. 
This  result  indicates  that  as  few  as  six  and  probably  even  four  angles 
would  be  sufficient  for  reasonable  accuracy. 
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Appendix  A 
PRADTN 


o  o  o  o  o 


AIBFiC  PRADTN  FULIST, DECK, REF  PRADCOOO 

SUBROUTINE  RADTN  PRAOOOK 

CMPILEO  OCTOBER  7,  1965  HBL 

PLANE  CHARACTERISTIC  TRACE  WITH  DOUBLE  GAUSSIAN  INTEGRATION  PRAD0030 

♦  ♦♦•♦♦♦♦♦♦♦♦♦♦♦♦♦A****************** *♦*****«>  ********** *******#***>*****PRADC04C 

*  SPUTTER  COMMON  **PRAD005C 

♦PRAD0060 

COMMON  LMOAI 37) ,  NR  ,  NSMLR  ,  IA  ,  IB  ,  ICA  ,  ICB  .PRA00070 

1  KMAX  ,  BLANK1,  BLANK2,  BLANK3,  IAP1  ,  IBP l  ,  ICAP1  ,  IC8P1  .PRAD0080 

2  II  ,  IG  ,  NRAD  ,  BLANK*,  IAM1  ,  IB Ml  ,  ICAMl  ,  1CBM1  ,PRAD0090 

3  II PI  ,  IGM1  ,  I ALPHA,  BLANK5,  TH  ,  TMAX  ,  8LANK6,  OELPRT , PRADO 1 00 

4  FREQ  ,  CNTMAX,  AR  ,  ASMLR  ,  PUSHA  ,  PUSHB  ,  BOILA  ,  BOILB  .PRADOllO 

5  CVA  ,  CVB  ,  SLUG  ,  ALPHA  ,  KVA  ,  HVB  ,  HCA  ,  HC8  ,PRAD0120 

6  EMINA  ,  EMINB  ,  CA  ,  CB  ,  GA  ,  GB  ,  GL  ,  GR  .PRAD0130 

7  RHOL  ,  RHOR  ,  EP10  ,  EPSI  ,  RIA  ,  RIB  ,  RDIA  ,  RDt B  .PRAD0140 

8  RPIA  ,  RPIB  ,  RPDIA  ,  RPDIB  ,  TPRINT,  TA  ,  TB  ,  TC  PRAD0150 

COMMON  TO  ,  TE  ,  DTH2  ,  DTH2P  ,  DTH1  ,  OTRMIN,  OTMAX  ,PRAD0160 

1  DTMAX1,  DTMAX2,  0TMAX3,  DTR  ,  SNITCH,  CO  ,  CHIN  ,  DELTA  ,PRAD0170 

2  GAMA  ,  UCRIT  ,  SIGMAQ,  AC  ,  ACQ3T4,  CNVRT  ,  SUMRA  ,  SUMRB  .PRAD0180 


LMOAI 37) • 
,  BLANK1, 
,  IG  , 
,  IGM1  , 
,  CNTMAX, 
,  CVB  , 
,  EMINB  , 
,  RHOR  , 
,  RPIB  , 
TO  , 


NR 

BLANK2 
NRAD 
I  ALPHA 

AR 

SLUG 

CA 

EPIO 

RPDIA 

TE 


DTMAX1, 
GAMA  , 


DTMAX2,  0TMAX3 
UCRIT  ,  SIGMAQ 


3 

ROI A 

,  ROIAMI, 

ROIB  , 

R0IBP1 

,  GMS 

.  SI 

,  S2 

• 

S3 

, PRAD0190 

4 

S4 

,  S5 

♦ 

S6  , 

S7 

,  SB 

,  S9 

,  sio  , 

Sll 

, PRAD0200 

5 

S12 

.  S13 

1 

S14  , 

S15 

,  S16 

,  S17 

,  S18  , 

S19 

, PRAD0210 

6 

S20 

.  EO 

1 

H>  t 

TAU 

,  ZERO 

,  R 

(152), 

0ELTARU52) 

»PRAD0220 

7 

ASQ 

(1521, 

RD 

(152) 

,  VO 

(152), 

ROD 

(152), 

SMLR 

(152) 

, PRAD0230 

8 

DELR 

I  37), 

P 

(152) 

,  PI 

(152), 

PB 

(152), 

PB1 

(152) 

PRA00240 

COMMON 

P2 

(152) 

,  SV 

(152), 

RHO 

(152), 

THETA 

(152) 

•  PRA002  50 

1 

U 

(152), 

E 

( 152) 

,  El 

(1521 , 

EK 

(152), 

A 

(152) 

, PRA00260 

2 

V 

(152), 

G 

(152) 

.  o 

(1521, 

C 

(152), 

X2 

(152) 

,PRAD0270 

3 

X3 

(152), 

X4 

(152) 

.  X5 

(152), 

X6 

(152), 

X7 

(152) 

.PRAD0280 

4 

SMLA 

(152), 

SMLB 

( 152) 

,  SMLC 

(152), 

SMLD 

(152), 

SMLE 

(152) 

, PRAD0290 

5 

EC 

(152), 

ER 

(152! 

,  SMLQ 

(152!  , 

SMLH 

(152!, 

B  (GA 

(152) 

•PRA00300 

6 

BIGB 

(152), 

CV 

(152) 

,  BC 

(152) , 

BR 

(152), 

CHIC 

(152) 

,PRAD0310 

7 

CHIR 

(152), 

CAPAC  (152) 

,  CAPAR 

(152), 

CRTC 

(152), 

CRTR 

(152) 

, PRA00320 

8 

CRT  PC 

( 152), 

GOER 

(152) 

,  FEU 

(152), 

CAR 

(152), 

OKLM 

(  37) 

PRAD0330 

COMMON 

TECM 

(  37) 

,  EKLM 

(  37), 

ELM 

(  37), 

FCLM 

(  37) 

,PRAD0340 

1 

FRLM 

(  37), 

HLM 

(  37) 

,  QLM 

(  37), 

AMASNOI  37), 

CHRNO 

(  37) 

, PRAD0350 

2 

2P1 

(  37), 

ZP2 

(  37) 

,  SOL  10 

(  37), 

ECHCK 

(  37), 

RK 

(104) 

, PRAD0360 

3 

4 

RL 

HEAD 

(  37), 

(  12), 

RHOK 

MAXL 

( 104) 

•  ROK 
,  MAXLM 

(104), 

THETAM  104) , 

TEMP 

(  16) 

.PRAD0370 

PRAD0380 

C*  **PRAD0390 

C  ****************************************************** *********«******«rRA00400 


DIMENSION  Q3(1),TG(1),H2( ll,Qlil),X8(l) ,SUMX3i 1 I *SUMX4( 1 I 
DIMENSION  H4(1),Y2(1)*H(1), SUMX2I 1),Q2(1) 

DIMENSION  037(1),  038(1),  H3fli 

COMMON  /LINOLY/  HNU, SGNL , IHNU, NHNU, HNUP »NT, IM , I N, DHNU, THICK, NY 
COMMON  /CNTRL/  SCYCLE,  JMULT 
COMMON  /DAVIS/  ICX,  ICY 

COMMON  /TQ/  Q INTI I  300) ,  QINT2I300),  TITLE! 12) 

EQUIVAL  ENCEI SML A»H4) , I SMLD,Y2) 

EQUIVALENCE  I BC,TG ), I  3 IGB , H), |CRTR,SUMX2» , (CHIC, SUHX3) 
EQUIVALENCE  ( SHLH, X8 ) , ( CAR ,Q37 ) , (CHIR, Q38) , ( SMLC ,H3) 
EQUIVALENCE  ( AC03T4, TRDBG ) , ( S 1 2, ED ITMF ) 

EQUIVALENCE  ( EC ,Q1 ) , (X7 ,H2) , ( S IGA, SUM a4» . (G0FR.Q3 ) 


PRA00410 

PRAD0420 

PRA00430 

PRAD0440 

PRAD0450 

PR4D0460 

PRA00470 

PRA00480 

PRAD049P 

PRAD0500 

PRAD0510 

PRAD0520 

PRAD0530 

PRAD0540 
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oooooooonooooor>nor>oooooooooooooooonooo 


t*****************;***********************«***********»****************pRA005;0 

♦PRAD0560 


EOITMf  SAME 

AS 

S12 

♦PRA00570 

H 

SAME 

AS 

BIGB 

•PRA00580 

AS  X7 

H3 

SAME 

AS 

SMLC 

PRA00590 

•PRA00600 

H4 

SAME 

AS 

SMLA 

•PRA00610 

01 

SAME 

AS 

EC 

♦PRAD0620 

03 

SAME 

AS 

GOFR 

•PRA00630 

037 

SAME 

AS 

CAR 

•PRA00640 

Q3S 

SAME 

AS 

CHIR 

•PRAD0650 

SUMX2 

SAME 

AS 

CRTR 

•PRA006&0 

SUMX3 

SAME 

AS 

CHIC 

•PRAD0670 

SUMX4 

SAME 

AS 

BtCA 

•PRAD0680 

TG 

SAME 

AS 

BC 

♦PRA00690 

TROBG 

SAME 

AS  AC03T4 

•PRA00700 

X8 

SAME 

AS 

SMLH 

•PRA00710 

V2 

SAME 

AS 

SMLO 

♦PRAD0720 

»**•*•< 

»•*** 

***************** 

•PRAD0730 

SUMRHO,  CSOO,  XSQO.  V« 

YSQD 

,  02*  NOT  USED 

PRA00750 

PRA00760 

PRA00770 

PRA007S0 

t**********************«************«<:»***«*****************PRAD0790 

PI  A 

N  E  S 

0 

N  L  V 

♦PRA00800 

•PRAD0810 

»*«***> 

***************** 

•PRA00820 

B  0  U  N  0  A 

R  V 

C  0 

N  0  I  T  I  0  N  S 

•PRAD0840 

PRA00850 

******* 

•PRA00860 

L.H.  BLACK BOOT 

IF  SOLIO  ON 

LEFT 

PRAD0880 

PRAD0890 

NO  VAPOR  ZONE  HAS  FLUX  OUT  FROM  SOLIO 

NTIMES-BOILB 
IM»l BM1 
IN*  I A 

IF<£Pl(26KEQ.O.)  CO  TO  15 

SAVE  STUFF  FROM  EIONX  FOR  NONEQ  AND  RESET  IN  OR  IM 
IF  (PUSHA  .IT.  O.OI  GO  TO  100 
IM  «  NR  -  1 


PRAD0900 

PRA00910 

PRA00920 

PRAD0930 

PRAD0940 

PRA00950 

PRA00960 


WSZ2 

WSZ3 

WSZ4 

HSZ5 


■BC(IM*1) 
•BRUNEI  > 

*  CRTCI IMUI 

•  RHO(  I M 1 ) 


100 


CO  TO  15 
IN  *  NR 

MSZ2  *  BC< IN-1) 
MSZ3  *  BR( IN-11 
WSZ4  *  CRTCI IN-1) 
WSZ5  *  RHO I  IN-1) 


’6 


<->  O  U  O  O  O  O  O  O  ooooooo  ooo 


15  CONTINUE 
IMP1«IM*1 
I NM1» I N-l 

CALL  CVCHK  (K000FX) 

IF  ( I HP1-I N I  190,190,125 

NO  VAPOR  20NES 

190  X2UMP1  )  *  1.0283612  *  MINI  *  < THETA! I N) **4  -  THETA!  IHPl )**4) 
ER(IH)«-X2(IMP1) 

GO  TO  1300 
125  IR*IN 

THT  AMX* .025 

IF  ( I ALPHA-li  130,150,130 
130  SI  *  13.0130 
CALL  UNCLE 
150  DO  ISO  I«IN,IM 
X3(I)*0. 

X5( 1 1*0. 

x5in*o. 

X6(t>»0. 

CRTR(I)*0. 

SET  UP  FOR  KAPPA  INTERPOLATION 

Ql(I)*THETA(n**4 
Q37I I )* ALOG! THETA! I ) ) 

Q38(  DIALOG!  SV<  III 

FIND  IK,  RIGHTMOST  ZONE  WITH  THETA  GREATER  THAN  0.05  EV 

IF  (THETA! I l-THTAMXI  160,160,150 
150  THTAMX*THETA(I) 

160  IF  (THETA! 1 1-0.051  i  80,170 
170  IP*  I 
180  CONTINUE 

IF  ( THT AMX  .LT.  THETA!  IB))  THTAMX  *  THETA! IB) 


PRA00990 

PRAD1000 

PRAD1010 

PRA01020 

PRA01030 

PRA01050 

PRAD1050 

PRAD1060 

PRA01070 

PRA01080 

PRA01090 

PRA01100 

PRAOillO 

PRA01120 

PRA01130 

PRAD1150 

PRAD1150 

PRAD1160 

PRA01170 

PRA01180 

PRA01190 

PRA01200 

PRA01210 

PRAD1220 

PRAD1230 

PRAD1250 

PHA01250 

PRA01260 

PRAO1270 

PRA01280 

PRAD1290 

PRAD1300 

PRA01310 

PRAD1320 

PRAD1330 

PRA01350 

PRA01350 

PRA01360 

PRA01370 


BEGIN  FREQUENCY  LOOP 


♦PRA01390 
•PRA01400 
*PR AD1410 


200  HNUP  *  3.E3  PRAD1430 

PRA01440 

SET  UP  MAX  FREQ  BOUNOARY  PRAD1450 

PRAD1560 

HNUP4  =  8.1E13  PRA01470 

I HNU*  1  PRA01480 

DO  210  l*  l  N*  IMP1  PRA01490 

210  SUMX2! I  1  =  0.0  PRAD1500 

IF  (KMAX.EQ.OI  GO  TO  280  PWAD151P 

C  PRA01520 

C  THIS  COOING  WONT  WORK  IF  HNU  NOT  EVALUATED  PRAD1530 


37 


OOO  UUO  OOO  O  O  <J  o 


c 

230  CALL  KAPPA! IN, INI 
HNU4»MNU**4 
OHNUP  *  DHNU 
OHNU  *  HNUP  -  HNU 

MERGE  GROUPS  WITH  HNU  MORE  THAN  CB  TIMES  LARGEST  THETA 

IF  I THTAMX-  HNU/CBI  240,300.300 

REJECT  TAPE  IF  MORE  THAN  HALF  OF  GROUPS  MERGE 

240  IF  (IHNU+IHNU-NHNUI  260,250,250 
250  IF  ! AMODICB, 1, I  .EQ.  0.51  GO  TO  260 
Sl-13.0250 
CALL  UNCLE 
260  00  270  I-IN.IM 
BETA*HNU/THETAIIt 
BETAP*HNUP/THETA!II 
OFBxPLNKUT ( BETA, BETAP ) 

IF  (OFB.EQ.O.I  GO  TO  270 
TEMP! 1 I *DF8*Q1I I ) 

EMBl-EXPi-BETAI 
EMB2*EXPI— BETAPI 

TEMP! 21 *DF8*0.0384974/Q1! I l*!HNU4/( 1.0-EMB1I 
1*EMB1-HNUP4/I 1.0-EMB2I*EMB2I 

FORM  NUMERATORS  ANO  DENOMINATORS  OF  MERGED  KAPPAS 

X6i I l*X6( I •♦TEMP! 1 ) 

X41 I l*X4! I ) ♦TEMP! 21 
X54  I  I«X5( 1 1 ♦CAP AC! I I*TEMP1 II 
X3I I )*X3( I J+TEMP! 21/CAPAR! I) 

270  CONTINUE 

IF  (GL  .IT.  I.  .OR.  I HNU  .EQ.  II  GO  TO  275 
MERGE  FREQUENCY-OEPENOENT  EXTERNAL  INPUT  INTENSITIES 
NMU  «  LMOA! 371 
IQNT  «  NMU  *  ! I HNU  -  21 
00  272  I  -  I,  NMU 

IQNT1  -  IQNT  ♦  I 
IQNT2  -  IQNTl  ♦  NMU 
IF! IHNU.GT .2 1  0HNUP*1. 

272  QINT1! IQNT2)  *  QINT11 IQNT2»*0HNU  ♦  QINTII IQNT1I*0HNUP 
275  HNUP-HNU 

IHNU* IHNU* 1 
HNUP4«HNU4 

IF  (THTAMX-  HNU/CBI  220,310,310 

MONOFREQUENCY  CALCULATION 

280  NMNU* 1 

CALL  KAPPA  ( IN, IM) 

DO  290  l=IN,IM 
X5! 1 1*1. 

290  X6(II>01( I) 


PRA01540 

PRADI550 

PRA01560 

9/29/65 

PRAD1570 

PRAD1580 

PRA01590 

PRA01600 

PRADI610 

PRA01620 

PRAD163Q 

PRA01640 

PRADI650 

PRA0I660 

PRA01670 

PRA01680 

PRA0I690 

PRA01700 

PRA017I0 

PRAD1720 

PRADI730 

PRA01740 

PRA0I750 

PRADI760 

PRA01770 

PRA01780 

PRADI790 

PRAD1800 

PRA0I810 

PRA01820 

PRA0I83O 

PR ADI  340 

PRA01850 

PRA0I860 

9/29/65 

PRA01880 

PRA01890 

9/29/65 

PRA0I9I0 

PRA01920 

PRAD1930 

9/29/65 

9/29/65 

PRA01950 

PRA01960 

PRADI970 

PRAD1980 

PRA01990 

PRA02000 

PRA02010 

PRA02020 

PRA02030 

PRA02040 

PRAD2050 

PRA02060 
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ooo  o  onn  o  ri  o 


DF8=1.0 
H!<U  *  .001 
ICX*IR 

IF  (GL  .GT.  0.01  ICX  *  IN 
ICY*IN 
GO  TO  480 

3C0  IF  IIHNU-1}  550.370.260 

FORM  MERGED  KAPPAS 

310  OC  350  I*IN, IM 

IF  ( X6(  1 1 )  320.350.330 
320  Sl~13.0320 
CALL  UNCLE 

330  CAPARI I l*X4! I )/X3( I ) 

340  CAPACi II*X5<II/X6i II 
350  CONTINUE 
HNUP  *  3.E3 
HNUP4  *  8.1E13 
OHNU  -  HNUP  -  HNU 
I HNU* IHNU-1 

IFIGL.LT. 1. »  GO  TO  370 
DO  355  1*1, NMU 
IQNT2  *  10NT  ♦  l  ♦  NNU 
355  QINTIIIQNT21  *  QINT1I  IQNT2I/DHNU 
GO  TO  370 

TYPICAL  GROUP  CALCULATION  OF  SOURCES 

360  CALL  KAPPA!  IN.IM) 

DHNU*HNUP-HNU 
HNU4*HNU**4 

370  IF  (GL-l.l  390.380,380 
380  IF  (HNU.NE.ROKI IHNU+52IJ  GO  TO  490 
IF  ( GL.NE. FLOAT ( NHNU I )  GO  TO  490 
ALCULATE  ICX,  ICY 
390  ICX*IN 
ICY  *  IN 

IF  <GL  .IE.  0.)  GO  TO  395 
ICX  «  IM 

00  392  I  *  IN,  IM 

OFB  *  PLNKUTIHNU  /  THETA! II,  HNUP  /  THETA! 1 1 1 
392  X6!  I  )*OFB*QH  1 1 
GO  TO  480 

395  00  470  I-IN.IR 

BET A* HNU/ THETA! I J 

AVOIO  CALCULATION  OF  OFB  LESS  THAN  IE-5 

IF  ! BET A-19.0  )  400,410,410 
400  8ETAP*HNUP/THETAI l I 
C  6MB2=EXP!-BETAP) 

IF  1BETAP-0.01)  410,410,460 
410  IF  (  ICX-IRI  430,420,420 
420  ICX*I-’ 


PRA02070 

PRA02080 

PRA02090 

PRA02100 

PRA02110 

PRAD2120 

PRA02130 

PRA02140 

PRA02150 

PRA02160 

PRA02170 

PRA02180 

PRA02190 

PRAD2200 

PRA02210 

PRA02220 

PRA02230 

PRAD2240 

PRAD2250 

PRAD2260 

9/29/65 

9/29/65 

9/29/65 

9/29/65 

PRA02270 

PRA02280 

PRA02290 

PRA02300 

PRA02310 

PRA02320 

PRA02330 

PRA02340 

PRA02350 

PRAD2360 

PRAC2370 

PRA02380 

PRAD2390 

PRA02400 

PRA02410 

PRA02420 

PRA02430 

PRA02440 

PRA02450 

PRA02460 

PRAD2470 

PRA02480 

PRA02490 

PRA02500 

PRA02510 

PRA02520 

PRA02530 

PRA02540 

PRA02550 

PRA0.7560 


non  r>  r>  n  non  o  n  n  <■>  o  o  r>  o  n  o  o  n  r>  o  o 


430  IF  II-ICYI  440,440,450 
440  ICYUCY+1 

ICX  IS  IN0EX  OF  LAST  ZONE  WITH  SIGNIFICANT  SOURCE 
ICY  IS  INDEX  OF  FIRST  ZONE  WITH  SIGNIFICANT  SOURCE 

450  X6( I ) *0  .0 
X  5( I ) *0.0 
GO  TO  470 

FORM  SOURCES  X6  AND  X5 

460  0FB*PLNKUT( BETA, BET API 
X6(!I«0FB*Q1U) 

TEMPI  2 >  *0. 0364974/011  I )*IHNU4/ 1  £  XP I BETA 1-1 . 0) 

1  -HNUP4/(1.0-EMB2I^ENB2I 
X5m«DIB»TEMP(2J 
ICX*IR 

470  CONTINUE 

480  IF  I INMli  490,520,500 

SET  6LACK800Y  CONDITION  FOR  !A  GREATER  THAN  1 

490  Sl=13.0490 
CALL  UNCLE 

500  OFB  *  PLNKUT  I HNU/THETAI INM1 » , HNUP/THETAI INM1 I » 

X6IINM1I  *  OFB  ♦  THETAIINM1 1**4 

SET  BLACKBOOY  CONDITION  IF  DESIRED  FOR  IMP! 

520  IF  (GL.NE.O.SI  GO  TO  530 

OFB  *  PLNKUTIHNU  /  THETAIIHPII,  HNUP  /  THETAIIMP1M 
X6IIMP1I  »  OFB  *  THETA! IMP  1 J**4 
530  031=0.0 

FORM  ROSSELANO  AND  PLANCK  OPTICAL  DEPTHS 

DO  590  1  =  1 N,  IM 
IF  (CAPARIIII  550,550,540 
540  IF  (CAPACUI)  550,550,560 
550  $1*13.0550 
CALL  UNCLE 

CHOOSE  ALL  ROSSELANO  IF  SOLIO  10  IS  POSITIVE 

560  IF  (SOLIOI 10I.EQ.0.I  GO  TO  570 
HI D-CAPARC ll/SV(II 
GO  TO  580 

570  Hm*CAPACm/SV(!J 
580  H2I  n=HU>*OELTAR(  I  J 

IF  (ALPHA  .GT.  l.»  GO  TO  586 
H  3 1  I >  *  CAPAR(I)  *  GIII 
GC  TO  588 

kVEAT.  ASYNCHROMSMS  IN  SV  ANO  OFLTAR  LEAD  TO  ERRONEOUS  FLUCTUATIONS 


PRAD2570 
PRAD2580 
PRAD2590 
PRAD2600 
PRAD2610 
PRAD2620 
PRAD26 SO 
PRAD2640 
PRAD2650 
PRAD2660 
PRAD2670 
PRAD2680 
PRAD2690 
PRAD2700 
PRAD2710 
PRAD2720 
PRAD2730 
PRAD2740 
PRAD2750 
PRAD2760 
PRAD2770 
PRAD2780 
PRAD2790 
PRAD2B00 
PRAD2810 
PRA02820 
PRAD2830 
PRAD2840 
PRAD2850 
PRAD2860 
PRAD2870 
PRAD2880 
PRAD2890 
PRA02900 
PRAD2910 
PRA02920 
PRA02930 
PRAD2940 
PRA02950 
PRAD2960 
PRAD29', 
PRAD2960 
PRAD2990 
PRA0300G 
PRA0301G 
PRA03020 
PRAD3030 
PRAD3040 
PRA03050 
PRAD3060 
PRAD3070 
PRA03080 
PRAD3090 


IN  H3.  THIS  CaN  BE  FIXED  BY  SUBSTITUTING  G  IN  PLANES,  BUT  SPHERE SPR AD 3 l 00 


WILL  STILL  HAVE  THIS  TROUBLE. 


PRAD3110 


U  O  U  OOO  o  o  u  o 


586  H3IIJ-CAPARI I)/SVII)*0£LTAR{I) 
588  031=Q31*H3I  l » 

Q3( I-*1)=Q31 
HIII<0.5*H(|| 

821  l)=0.5*H2IIi 
H3ll)=0.5*H3II) 


ZERO  D I FFUSU  IN  INDICATORS  AND  X2 


X  2 1  I  >*0.0 
X3 1  I >  *0  .0 
X4Il)=G.O 
590  R8C(I)=0.0 
X2I IMP1 )=0.0 
X3I IMP1 )=0.0 
X4 1 I MP1 )=0 .0 
IF  ( icr  .GT. 


GO  TO  990 


STEP-LINEAR  CRITERION  AT  ICY 
UNCONDITIONAL  STEP  AS  BOUNDARY  CONDITION  IF  ICY  *  IN 

IF  (  ICY-IN)  675,600, 610 
600  Y2( IN)=X6( IN) 

TGI  I N )  =  0.0 
GO  TO  620 

610  TEMPI 1)=H3( ICY-1 ) ♦H3I ICY ) 

TGI ICY)=X6CICY)/TEMP(1) 

Y2I ICY) =TG! ICYI»H3< ICY-i) 

FORM  Y2  AND  TG  SET  X3=-l  IF  A  DIFFUSION  CRITERION  MET  USING  HCB 


620  ICXM1*ICX-1 

IF  I  ICY  .GT.  ICXM1)  GO  TO  67 2 
DO  670  I=ICY« ICXMI 
TEMP{1)=H3I  I  +  D+H3I  I) 

IF  I ANAX1IX6I I),  X6II+1))  .LE .  0.)  GO  TO  640 
IF  (AMINIIH3I I).  H3II+1))  .GT.  AC)  GO  TO  650 
IF  IABSIIH3I  I ) — H 31  I ♦ 1) I/TEMPI  1)  )  .GT.  GA)  GO  TO  640 
IF  I ABSI  (  X6I  I  )-X6<  1 ))/{  X6I I )  *X6»  I *1) )  )-GA)  650,650,640 
640  TGI I+1)=0.0 
GO  TO  670 

650  TGI  I ♦11  =  1X61  I*1)-X6I  I))/TEMPU) 

TGI !♦ I )  =  1 0 II I ♦  ll-Ql I 1 ) ) /TEMPI 1)*IX5(I>1)+X5 11)1/2.0 

Y2I I  ♦  1)  =  I  X6I  I  ♦  I » *H3I  I J+X6I  I)*H3t  I ♦  1 ) ) /TEMPI  1 ) 

IF  (ABSITGI  I*l))-HCB*Y2( I»l)l  660,670,670 
669  X3I I  M)  =-l  .0 
6'.  0  CONTINUE 

RADIATION  BOUNDARY  CONDITION  AT  ICX 
I  VACUUM  IF  ICX  «  IM  AND  GL  NOT  1/2) 

672  IF  UCX-IM)  680,690,675 
675  SI  =  13.0675 


PRAD31 20 

PRA03130 

PRAD3140 

PRAD3I5 

PRAD316 3 

PRAD3170 

PRA03180 

PRAD3190 

PRAO3200 

PRA03210 

PRA03220 

PRA03230 

PRA03240 

PRAD3250 

PRA03260 

PRA03270 

PRAD3280 

PRAD3290 

PRAD3300 

PRAD3310 

PRAD3320 

PRA03330 

PRA03340 

PRAD3350 

PRAD3360 

PRAD3370 

PRAD3380 

PRAD3390 

PRAD3400 

PRAD34 10 

PRAD3420 

PRA03430 

PRAD3440 

PRAD3450 

PRAD3460 

PRAD3470 

PRAD3490 

PRA03500 

PRAD3510 

PRAD3520 

PRA03530 

PRAD3540 

PRA03550 

PRA03560 

PRA03570 

PRAD3580 

PRA03590 

PRA03600 

PRAD3610 

PK„03620 

PRA03630 

PRA03640 

PRAD3650 

PRAD3660 


- ...  ..„,e  % 


o  o  o  o  o  o  n  o  o  o 


CALL  UNCLE 

PRAD3670 

680 

TEMPI  1 )  *H3( ICX+1) +H3I ICX) 

PRA03680 

TGI ICXfll«-X6(ICXt /TEMPI  11 

PRAD3690 

Y2I  ICX+1  J  —  TGI  ICX+1 )  *H3(  ICX+l) 

PRA03700 

GC  TO  700 

PRA03710 

690 

IF  (GL  .EQ.  0.5)  GO  TO  700 

PRA03720 

V2IIMP1MX6IICX) 

PRA03730 

TGI IMP1 )*0 .0 

PRAD3760 

PRAD3750 

EXTEND  TRANSPORT  REGION  BOUNDARIES.  IF  NEEDED.  TO  PROVIDE  HVft  MEAN 

PRAD3760 

FREE  PATHS 

PRAD3770 

PRAD3780 

700 

I*  IN*  1 

PRAD3790 

710 

IF  (X3II))  720.760,730 

PRAD3800 

720 

1*1  +  1 

PRAD38 10 

IF  II-ICX-1)  710,760,820 

PRAD3820 

730 

Sl*13.0730 

PRA03830 

CALL  UNCLE 

PRAD3860 

760 

J*I-1 

PRAD3850 

750 

IF  (031 l  )-Q3( JJ-HVB)  760,760,770 

PRAD3860 

760 

X6IJ)*-1.0 

PRAD3870 

J*J-1 

>RAD3880 

IF  (J-IN)  770,750,750 

PRAD3890 

770 

1*1  +  1 

PRAD3900 

IF  II-ICX-1)  780,780,820 

PRAD3910 

780 

IF.'  .1X31  II)  790,770,730 

PRAD3920 

790 

J*l 

PRAD3930 

800 

IF  (031 J 1-031 I -1 l-HVB )  810,810,720 

PRAD3960 

810 

X6(J)*-1.0 

? RAD3950 

J*J+l 

PRAD3960 

IF  IJ-iCX-1)  800,720,720 

PRAD3970 

820 

I* I N+  1 

PRAD3980 

PRA03990 

TEST  TO  FORM  TT  ANSPOR T  REGIONS 

PRA06010 

830 

I  AX*  I N 

PRAD6020 

860 

IF  (X3IH)  850,860,730 

PRAD6030 

850 

IF  1X6(  11)  860.87C. 7 1 1 

PRAD6060 

PRAD6050 

REMOVE  ONE  ZONE  DIFFUSION  REGION 

PRAD6070 

860 

1*1  +  1 

PRAD6080 

IF  I I-ICr  -1)  860,950 , 950 

PRAD6090 

870 

1  =  1  +  1 

PRAD6100 

IF  :  l - 1 CX- 1 1  890,950,950 

PRAD6110 

880 

It  (X 3(H)  890,860,730 

PRAD6120 

89C 

».  X6I I ) )  860,900,730 

PRAD6130 

;f>0 

fp  *1-3 

PRAD6160 

GO  TO  960 

PRAD61 50 

910 

IF  I X3I I ) )  920,960,730 

PRAD6160 

920 

IF  (X6(U‘  960,930,730 

FORM  X2  FOR  DIFFUSION  ZONES  IN  ORDER 

PRA061 70 

930 

X 2 C I )  *  -1.37E12  ♦  TGI  I ) 

PRAD6190 

1*1  +  1 

PRAR620C 

IF  (I  (CX-1)  910,980,980 

PRAD62 10 

42 


o  o  o  o  o  o  uo  o  u  o  o  o  u  o  o 


PRAD4220 

PRAD4240 
PRAD4250 
PRA04260 
PRAD4270 
PRAD4280 
PRAD4290 
PR AD4300 
PRAD43I0 


GO  TO  950  PRAD4330 

PRA04340 

OPTIONAL  EOIT  OF  X2  ETC.  PRA04350 

PRAD4360 

990  IF  ( EOITMFI  1020,1020,1000  PRAD4370 

1000  IARG1*SQL 10(18) +0.001  PRA043 80 

IARG2*E0ITMF+0.0Q1  PRA04390 

IF  (MOD( I ARG1 , I ARG2) )  1020,1010,1020  PRA04400 

1010  WRITE  (3)  HNU ,  IN,  IN,  IMP1,  SOLIO(IS),  TH,  DHNU  PRA04410 

IN  SPHERICAL  VERSION,  REOIT  IS  GIVEN  RHO.  PRAD4420 

REPLACEO  HERE  BY  CAPAR.  PRA04430 

WRITE  (3)  <CU),  I  =  IN,  IMP1 )  »  (H3(I),  l*IN,IM),  (X6(I),  I *1 N, I MP1 ) .PRAD4440 
1(X2(  I ),I  =  IN,  IMP1),  PRA04450 

2  (CAPAR(I),  I*IN,lMPl),  (THETA(I),  1*IN,IMP1),  (EKI),  I*IN,IM)  PRAD4460 
XX*~2.0  PRAD4470 

WRITE  (31  XX,XX,XX,XX,XX,XX,XX  PRA04480 

BACKSPACE  3  PRA04490 

JNUL  T*  1  PRAD4500 

1020  OC  1030  I=JN.IMP1  PRA04510 

SUMX2( I )*SUMX2( I )  +  X2(  I)  PRAD4520 

1030  CONTINUE  PRA04530 

PRA04540 

ADVANCE  FREO,  STORE  EMERGENT  FLUX,  TEST  FOR  COMPLETION  OF  GROUPS  PRAU4550 

PRA04560 

HNUP=HNU  PRA04570 

HNUP4=HNU4  PRAD4580 

I HNU  *  IHNU  ♦  1  PRAD4590 

IF  (IHNU-NHNU)  1040,1040,1060  PRAD4600 

1040  CALL  OVCHK  (KOOOFX)  PRA04610 

GO  TO  (1050,360),  KOOOFX  PRA04620 


DO  TRANSPORT  TO  1M  IN  REGION  OF  NO  SOURCE 

940  I  AX* I 

GO  TO  860 
950  18X*IM 

960  CALL  TRANS( I  AX, IBX) 

IF  (IBX-IM)  970,990,990 
970  I  *  I BX  +  2 
GC  TC  9  30 

980  IF  ( I  .GT.  IM)  GO  TO  990 
I  AX*  I 
GO  TO  950 

OPTIONAL  EOIT  OF  X2  ETC. 

990  IF  ( EOITMF)  1020,1020,1000 
1000  IARG1*SQL 10(18) +0.001 
IARG2*EOITMF+0.001 

IF  (MOD( I ARG1 , I ARG2) )  1020,1010,1020 
1010  WRITE  (3)  HNU.  IN,  IM,  IMP1,  SOLIO(IS),  TH,  OHNU 
IN  SPHERICAL  VERSION,  REOIT  IS  GIVEN  RHO. 

REPLACEO  HERE  BY  CAPAR. 


HNUP=HNU  PRAD4570 

HNUP4=HNU4  PRAD4580 

IHNU  *  IHNU  ♦  1  PRA04590 

IF  (IHNU-NHNU)  1040,1040,1060  PRAD4600 

1040  CALL  OVCHK  (KOOOFX)  PRA04610 

GO  TO  (1050,360),  KOOOFX  PRA04620 

,♦*»**♦,♦*»**»»*»*»«.**♦********♦***»**,**♦***»»*♦♦****»**♦*++♦****» *++,PRAD4t 30 

♦PRAD4640 

ENO  FREQUENCY  LOOP  +PRAD4650 

♦PRA04660 

******* ******************** ************ ***** **********  ******  ******* ****PRA046 70 

1050  Si  =  13.1050  PRAD4680 

CALL  UNCLE  PRAD4690 

1060  SUMX2 ( l NMl)  »  0.0  PRAD4700 

00  1070  l*INHl, (MP1  PRA04710 

X2( I )  =  SUMX2 ( I )  PRA04720 

1070  ER (  I  )  *  SUMX2 (  I )  -  SUMX2U  +  1)  PRA04730 

C  PRA04740 

C  FORM  MGNOFREOUENCY  QUANTITIES  ANO  FIND  MIN  TIME  STEP  PRAO4750 


43 


o  u  o  u  o  o  o  o  uut> 


c 

MSB  -  0.0 

00  1075  I  »  1.  MAXLM 

1075  MSB  -  MSB  ♦  ELMI  1 1 
DTR1-1.E10 
DTR2»i.E10 

IF  ( KMAX. EQ.O)  GO  TQ  1080 
CALL  KAPPA! IN. IN) 

1060  00  1230  I *IN, IN 

IF  ROSS  IS  ZERO  EXIT 


1090 


1110 

1120 

1130 

1140 


1160 


1165 

1170 


IF  (CAPAR(l))  1090.1090.1100 

Sl-13.1090 

oALL  UNCLE 

TEMPI 3)«CAPAR( I) 


IF  I SCLIOI 10) )  1110.1120.1110 

tempi  d>capar;:j 

GO  TO  1130 

TEMPI  1) -SQRTICAPARI I )*CAPAC( I ) ) 

IF  (TEMPI  1))  1090,1090,  1KC 
HI  I )«0. 5*TEMPI 1 )/SV( I ) 

H3II)  *  HII)  *  DEL TAR 1 1 ) 

IF  I .00 l-THETAI I ) )  1160,1230,1230 

IF  (H3I  n.GT.0.1)  GO  TO  1170 

IF  (ERI  D.EO.O.)  GO  TO  1170 

WSSB  -  Eli)  ♦  GUI 

IF  (TCLN( 37 )  .EC.  0.0)  GO  TO  1170 

IF  I MSBB  -  TELMI 37 )  *  MSB)  1170,  1165,  1165 

TEMPI 2)*.5«CV(I)*THETA( I )*G( I  )/ABS(ER( I )) 

GO  TO  1180 

TEMP<2)=»<.5U.5*H3(l)**2)*CVI  1)/<4.1132E12*TENP(3)*THETA(I>4*3) 
TEMPI 2)»TEMP(2)*TELM( 25) 

.  ***** ****** **************************************  ************** 


PR ADA  760 
PRAD4770 
PRAD4780 
PRA04790 
PRA04800 
PRAD48 10 
PRA04820 
PRAD4830 
PRA04840 
PRAD4850 
PRAD4860 
PRAD4870 
PRAD4880 
PRAD4890 
PRAD4900 
PRAD4910 
PRAD4920 
PRAD4930 
PRA04940 
PRAD4950 
PRAD4960 
PRAD4970 
PRA04980 
PRAD4990 
PRAD5000 
PRAD5010 
PRAD5020 
PRAD5030 
PRAD5040 
PRAD5050 
PRAD5060 
PRA05070 
PRAD5080 


•PRAD5100 

FIND  MINIMUM  TIME  STEP  PPRAD5110 

♦PRAD5120 

ft*********************** ******************* •*******«<>***«»* *****«**4**PRAD5 130 


1180 

IF  I  TEMPI  2) )  1230,1230,1190 

PRA05140 

1190 

CONTINUE 

PRAD5150 

IF  I  TEMPI  2 )— OTR 1 )  1200,1210,1210 

PRAD5160 

1200 

0TR2-DTR1 

PRAD5170 

IMN2-IMN1 

PRAD51 80 

OTR l»  TEMPI  2 ) 

PRAD5190 

I MNl» I 

PRAD5200 

GO  TO  1230 

PRA05210 

1210 

IF  I  TEMPI 2I-0TR2)  1220,  1230,1230 

PRAD5220 

1220 

OtR2-TbMPI2) 

PRAD5230 

IMN2-I 

PRA05240 

1230 

CONTINUE 

PRAD5250 

DTRMIN“0TRl 

PRAD5260 

E0UMN1 

PRAD5270 

PRA05280 

PRINT  MINIMUM  TIME  STEPS  BEiMEEN  EDITS 

PRAD5290 

PRAD5300 

44 


ooo  ooooo 


IF  I0TR1-TELMI26) )  1240, 1250t 1250 
1240  TELMI26)=DTR1 
TELH(27)*IMN1 
T  ELMI 28 )*DTR2 
TELMI  29IMMN2 
TELMI 30  >= S 0C I D( 18I  +  1.0 
1250  CONTINUE 

OETERHINE  IF  RADIATION  OR  HY0R0  WILL  SUBCYCLE 

IF  (DTRMIN-OTR)  1280,1300,1260 
1260  BLANK3*TH«-AMIN1(0TRMIN,GR*DTH2) 

IF  (S17)  1300,1270,1300 
1270  S9  «  1,0 
GC  TO  1300 

*#«4A*',«**lM***«*««******«******«*t**«*****************«***#**t********pRAD5460 

•PRAD5470 

REDUCE  TIME  STEP  *PRAD5480 

♦PRA05490 

♦♦♦♦♦♦•♦♦♦♦♦♦***********#4***#**»4*******************»***********#*****pRA05500 


PRA05310 

PRAD5320 

PRAD5330 

PRA05340 

PRA05350 

PRAD5360 

PRA05370 

PRAD5380 

PRA05390 

PRAD5400 

PRA05410 

PRA05420 

PRAD5430 

PRA05440 

PRA05450 


1280  NRAD*ZP 1 (18)/0TRMIN+1,0 
DTR*Z  PI I 18I/FLOATINRAO) 

IF  (NRAD-NT IMESI  1300,1300,1290 
1290  SU13.1290 
CALL  UNCLE 

ZERO  OUT  STRAY  QUANTITIES  FROM  PREVIOUS  CYCLES 
1300  DC  1310  I  *  IMP  1 ,  IG 
CAPARI 11*0. 

CAPACII)  -  0. 

X2IU1)  «  0. 

X3 1  I  ♦  1 1  =  0. 

X4( I ♦ 1)  »  0. 

SUHX2IIMI  -  0. 

SUMX3 II ♦ 1 I  «  0. 

SUMX4I 1*11  »  0. 

1310  ER( I »ll  »  0. 

IF  (ZP1I26S  ,EQ.  0.0)  GO  TO  1400 
:  RESTORE  GOOOIES  FOR  NONEQ 

.  (PUSHA  ,LT.  0.0)  GO  TO  1350 

BC  IMP1  )«WSZ2 

BRI IMP1 )=WSZ3 

CPTRI IMP1)  «  WSZ4 

RuOI  IMP1)»WSZ5 

GO  TO  1400 

•'  .3  BC(INMl)  «  WSZ2 
BRI INM1 I  =  WSZ3 
CRTRI  I N Ml I  «  WSZ4 
RHCI INM 1 )  *  WSZ5 
1400  RETURN 
END 


PRA05510 
PRA05520 
PRA05530 
PRA05540 
PRA05550 
PRA05560 
PRAD5570 
PRAD5580 
PRAD5590 
PRA05600 
PP.AD56 10 
PRAD5620 
PRA05630 
PRA05640 
PRA05650 
PRA05660 


PRA05680 


4< 


This  page  intentionally  left  blank. 
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UBFTC  PTRANS  FUL 1ST ,OECK,REF  PTRAOOOO 

SUBROUTINE  TRANSIN, M)  PTRAOOIO 

CCMPILEO  JULY  1,  1965  W8L  PTRA0020 

PLANE  CHARACTERISTIC  TRACE  WITH  DOUBLE  GAUSSIAN  INTEGRATION  PTRA0030 

DECK  OBLGAU  REQU1REO  FOR  INTEGRATION  COEFFS.  PTRA0040 

************************************************************ ***********PTRAOO 50 


SPUTTER  CONHON 


COMMON 
KM  AX 
II 

1 1  PI 
FREQ 
CVA 
EMINA 
RHCL 
RP I A 


1 

2 

3 

4 

5 

6 

7 

8 

C 

1 

2 


COMMON 


0TMAX1, 
GAMA  , 


LMDA1 37) , 
•  BLANK l, 
,  IG  , 
,  IGM1  , 
,  CNTMAX, 
,  CVB  , 
,  FMIN8  , 
,  RHOR  , 
,  RPIB  , 
TO  , 
0TMAX2, 
WCRIT  , 


NR 
BLANK2 
NRAO 
I  ALPHA 

AR 

SLUG 
CA 
EP  10 
RPOIA 
TE 

0TMAX3 

SIGMAQ 


NSMLR  , 
BLANK3, 
BLANK*, 
BLANK5, 
ASMLR  , 
ALPHA  , 
CB  , 
EPSI  , 
RPOIB  , 
DTH2  , 
OTR  , 
AC  • 


IA  , 
I  API  , 
I AN1  • 
TH  , 
PUSHA  , 
HVA  , 
GA  , 
R1A  , 
TPRINT, 
0TH2P  , 
SWITCH, 
AC03T4, 


IB 

IBP1 

IBM1 

TMAX 

PUSHB 

HVB 

GB 

RIB 

TA 

DTH1 

CO 

CNVRT 


ICA  , 
ICAP1  , 
!€AHl  , 
BLANK6, 
BOILA  , 
HCA  , 
GL  , 
ROIA  , 
TB  , 
DTKMIN, 
CMIN  , 
SUMRA  , 


♦♦PTRA0060 
♦PTRA0070 
ICB  .PTRA0080 
ICBP1  , PTRA0090 
ICBM1  , PTRAOIOO 
OELPRT . PTRAOl 10 
BOILS  *PTRA0120 
HCB  »PTRA01 30 

GR  , P/RA0140 

ROI 8  .PTRA0150 

TC  PTRA0160 

OTMAX  ,PTRAOl 70 
OELTA  ,PTRA0180 
SUMRB  .PTRA0190 


3 

ROIA 

•  R0IAM1, 

RO IB  ,  ROtBPl, 

GMS 

,  SI 

,  S2 

* 

S3 

,kTRA0200 

4 

S4 

.  S5 

, 

S6  ,  S7  , 

SB 

,  S9 

,  S10  , 

Sll 

,PTRA02i0 

5 

S12 

,  S13 

, 

S14  ,  S15  , 

S16 

,  S17 

.  SI 

8  , 

S19 

.PTRA0220 

6 

S20 

,  EO 

, 

FO  ,  TAU  , 

ZERO 

,  R 

(152), 

DELTAR( 152 ) 

»PTRA0230 

7 

ASQ 

( 1521, 

RO 

(152),  VO 

(152), 

RDD 

(152), 

SMLR 

(152) 

.PTRA0240 

8 

DELR 

(  371, 

P 

(152),  PI 

(152), 

PB 

(152), 

PB1 

(152! 

PTRA0250 

COMMON 

P2 

1152),  SV 

(152), 

RHO 

(152), 

THETA 

(152) 

.PTRA0260 

1 

W 

(152), 

E 

(152),  El 

(152)  , 

EK 

(152) , 

A 

(152) 

, PTRA0270 

2 

V 

( 152), 

G 

(152),  0 

(152), 

c 

(152), 

X2 

(152) 

,PYRA0280 

3 

X3 

(152), 

X4 

(152),  X5 

(152) , 

X6 

(152), 

X7 

(152) 

.PTRA0290 

4 

SMLA 

( 152), 

SMLB 

(152),  SMLC 

( 152) , 

SMLO 

(152)  , 

SMLE 

(152) 

*PTRA0300 

5 

EC 

< 152), 

ER 

(152),  SMLQ 

(152)  , 

SMLH 

(152) , 

BIGA 

(1521 

.PTRA0310 

6 

BIGB 

( 152)  , 

CV 

(152),  BC 

(152), 

BR 

(152), 

CHIC 

(152) 

.PTRA0320 

7 

CHIR 

( 152) , 

CAPAC  (  15  2,',  CAPAR 

(152) , 

CRTC 

(152), 

CRTR 

(152) 

.PTRA0330 

8 

CRTPC 

(152), 

GOFR 

(152),  FEW 

(152), 

CAR 

( 1 52  >  , 

OKLM 

(  37) 

PTRA0340 

COMMON 

TELM 

(  30,  EKLM 

(  37), 

ELM 

(  37), 

FCLM 

(  37) 

,PTRA0^50 

1 

FRLM 

(  37), 

WLM 

(  37),  QLM 

(  37)  , 

AMASNOt  37), 

CHRNO 

(  37) 

.PTRA0360 

2 

ZP1 

<  37), 

ZP2 

(  37),  SOLID 

(  37), 

ECHCK 

(  37), 

RK 

(104) 

.PTRA0370 

3 

RL 

(  37), 

RHOK 

(104),  ROK 

(10t) , 

THE  TAM  104)  , 

TEMP 

(  16) 

.PTRA0380 

4 

HE  AO 

(  12), 

MAXL 

,  MAXLM 

PTRA0390 

C* 


**PTRA0400 


C *************** ********************************.**♦#*** ****** ***********PTRA04 10 

DIMENSION  Q 3 ( 1 ) , T G < 1),H2( 1) ,01111 , X8( 1) , SUMX3II)  ,SUMX4C U  PTRA0420 

DIMENSION  H4(1),Y2I1), HID, SUHX2U), 02(11  PTRA0430 

D I  he  NS  ION  Q37ID,  Q38ID,  H3I  1 )  PTRA0440 

COMMON  /LINOLY/  HNU, SGNL , IHNU, NHNU.HNUP ,NT , I M , I N, DHNU, TH ICK ,NY  PTRA0450 

COMMON  /CNTRL/  SCYCLE,  JMULT  PTRA0460 

COMMON  /DAV IS/  ICX,  ICY  P7RA0470 

COMMON  /TO/  Q  ’NT  1 ( 300 ) ,  QINT2I300),  TITLEI12I  PTRA0480 

EQUIVALENCE  (SMLA,H4),(SHL0,Y2)  PTRA0490 

EQUIVALENCE  I BC, TGI , I  8  I GB , H) , I CRTR , SUHX2 ) , (CH IC , SUMX3)  PTRA0500 

EQUIVALENCE  I SMLH ,X 8 1 , I C AR, Q3 7 ) , ( CHIR, Q38 ) , I SMLC ,H3»  PTRA0510 

EQUIVALENCE  I AC03 T4 , TROBG ) , I S 1 2, EOITMF )  PTRAG520 

EQUIVALENCE  (EC, 011,1X7  , H2 ) , ( BIG A , SUMX4) , (GOFR , C3 )  PTRA0530 

L ******* ****************** ***************************** ****** ******** ***PTRA0540 


a 


f 


01  MENS  I CN  RR(40)  PTRA0550 

DATA  RR/2.113248E-01,7.886752E-01,1.056624E-01,3.94337&E-01,  PTRA0560 

1  1. 127017E-01, 5 .000000E-01 , 8. 872983E-01 ,3 .130600E-02  »  PTRA0S70 

2  2.222222E-01»2,464718E-01» 6.943180E-02 ,3.300095E-01»  PTRA0580 

3  6.699905E-01.9, 305482E-01, 1 « 20761 0E-02  »1 *07607 1 E-01 1  PTRA0590 

4  2. 184655E-01, 1.618513E-01 ,4. 691010E-02 t2 « 30765 3E-01,  PTRA0600 

5  5.000000E-01,7.692347E-01,9.530899E-01,5.557100E-03,  PTRA0610 

6  5, 522540E-02, 1.422222E-01 • 1. 840889E-0l» 1 .129063E-01 »  PTRA0620 

7  3#  37652QE-02, 1.693953E-01 « 3. 806903E-01 ,6 .193096E-01 ,  PTRA0630 

8  8.306047E-01,9.662348E-01,2.892400E-C3,3.055570E-02,  PTRA0640 


c 

9 

8.906520E-02, 1 • 44891 8E* 

■01, 

1 ,4982 5 IE-01 ,8*2 76980E-02/ 

PTRA0650 

♦PTRA0660 

c 

EOITMF 

SAME 

AS 

S12 

♦PTRA0670 

c 

H 

SANE 

AS 

BIGS 

•PTRA0680 

c 

H2 

SAME 

AS 

X7 

PTRA0690 

c 

H3 

SAME 

AS 

SMLC 

•PTRA0700 

c 

H4 

SAME 

AS 

SMLA 

♦PTRA0710 

c 

Q1 

SAME 

AS 

EC 

♦PTRA0720 

c 

Q3 

SAME 

AS 

GOFR 

•PTRA0730 

•+ 

037 

SAME 

AS 

CAR 

•PTRA0740 

c 

Q38 

SAME 

AS 

CHIR 

♦PTRA0750 

c 

SUMX2 

SAME 

AS 

CRTR 

♦PTRA0760 

c 

SUNK  3 

SAME 

AS 

CHIC 

♦PTRA0770 

c 

SUMX4 

SAME 

AS 

BIGA 

•PTRA0780 

c 

XG 

SAME 

AS 

BC 

♦PTRA0790 

c 

1R0BG 

SAME 

AS 

AC03T4 

♦PTRA0800 

c 

X8 

SAME 

AS 

SMLH 

♦PTRA0310 

c 

Y2 

SAME 

AS 

SMLO 

♦PTRA0820 

FEX,  FM,  SUMRHO,  CSQO, 


•PTRA0830 

,**»«9*«*««$««*««»*««***************#*******PTRA0S48 

PTRA0850 

PTRA0860 

XSQO,  Y,  YSQO,  Q2,  NOT  USED  PTRA0870 

PTRA0880 


C 

C* 

C 
C 
C 

c 
c 

C******4AA***#********#A******AA**A***********4P***********I*******A*****PTRA0900 


c 

c 

c 


PLANES  ONLY 


♦PTRA0910 

♦PTRA0920 

•PTRA0930 


C ********************************************************************** *PTRA0940 


C 

C 

C 


IAX«N 

IBX-M 

I  MIA 

INMl-lN-l 

IMP1  «  (M  ♦  1 

CALL  DVCHK ( KOOOFX ) 

GO  TO  (100,1101,  KOOOFX 
100  Sl*14.0100 
CALL  UNCLE 

no  i8xpi«iex*i 

1  ALPHA* ALPHA 

ERROR  IF  NOT  PLANE 
GO  TO  (130,120.120),  (ALPHA 


PTRA0950 

PTRA0960 

PTRA0970 

PTRA0980 

PTRA0990 

PTRA1000 

PTRA1010 

PTRA1020 

PTRA1030 

PTRA1040 

PTRA1050 

PTRA1060 

PTRA1070 

PTPA1080 

PTRA1090 
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120  SI* 14.0120 
CALL  UNCLE 

130  NY  *  LM0AI37)  -  1 

NNU  *  1NY  -  1)  *  (NY  ♦  2)  1 

NGS  *  NHU  ♦  NY  ♦  1 
JJ  *  0 

00  POSITIVE  ANGLES  FIRST 

140  1*1  AX 
F2*0.0 

IF  IAX*IN  TRANSFER  TO  150  TO  SET  SPECIAL  BOUNOARY  CONDITIONS 
IF  (IAX-IN)  360,150.180 

CALCULATE  BOUNOARY  SOURCE  INTENSITY 
150  IF  ( INN1 )  160,310,170 

SET  BLACKBOOY  CONDITION  FOR  PUSHER 

160  SW4.0160 
CALL  UNCLE 
170  F2*X6< INM1I 

GO  TO  310 

OIFFUS ION  BOUNOARY  CONOITIGN  AT  IAX 

180  IF  ( TGI  I )  .EQ.  0.1  Y2(  !)  *  X6( I-li 

TEHP|1I*H2( I- 1 )/RR (NNU I 
X95 I»*TG(l)*RRINHU) 

IF  (TEMPI 1)-1.E~2I  190,190,200 
190  F2X  *  Y2I I  —  1 1  -  TGII-ll  ♦  RRINMU) 

F2  *  IIV2II)  ♦  Y21 I-l ) )  *  0.5  *■  X6II-1)  -  F2X-F2X)  *  TEMPIl) 
GO  TO  250 

200  H4I l  — 1 )  *FR EX PI  — TEMPI  111 

F2*Y2( I )-X8( D+IX3I i l-RRI NNU) *TGI l-l»)*H4(I-l» 

GO  TO  250 

210  IF  ( TGI !-l )  .EQ.  0.)  Y2II-1I  *  X6II-1I 
X8I 1-1 ) *TG< I — 1 ) *RR I NMUI 

REGULAR  INTEGRATION  STEP  FOR  F2,  POSITIVE  MU 

220  IF  I  TGI  I )  .EQ.  0.)  Y2(I>  *  X6II-1) 

X8I I )*TG( I I*RR(NMU» 

TENPtl)«H2II-l)/RR(NMUJ 
IF  (TEMPI  1 1  —  l.E— 21  230,230,240 

230  F2  *  l(Y2(!I  +  Y2U-1I)  *  0.5  ♦  X6(I-1>  -  F2  -  F2>  *  TEMPI1J 
GO  TO  250 

240  H4I I- l I “FREXPI -TEMP ( 1 ) ) 

F2*Y2(  I  1 —  X  8  (  I)MIF?-Y2(  I  -  1 » ♦  X8<  l- 1 )  I  *H4 1 1- 1 » 

1*X8(I)-X8< I-l)  1-1) 

250  IF  (F2. IT. 0.0)  GO  TO  280 
260  SUMX3 ( I )*F2 


PTRA1100 
PTRA1110 
PTRA1120 
PTRA1130 
PTRA1140 
PTRA1150 
PTRA1 160 
PTRA1170 
PTRA1180 
PTRA1190 
PTRA1200 
PTRA1210 
PTRA1220 
PTRA1230 
PTRA1240 
PTRA1250 
PTRA1260 
PTRA1270 
PTRA1280 
PTRA1290 
PTRA1300 
PTRA1310 
PTRA1320 
PTRA1330 
PTRA1340 
PTRA1350 
PTRA1360 
PTRA1370 
P l RA1380 
PTRA1390 
PTRA1400 
PTRA1410 
PTRA1420 
PTRA1430 
F2XPTRA1440 
PTRA1450 
PTRA14  >0 
PTRA14 70 
PTRA14 JO 
PTRA1*  90 
PTRAlf  JO 
PTRAir 10 
PTRAl',20 
PTRA1530 
PTRA1540 
PTRA1550 
PTRA1560 
PTRA1570 
PTRA1580 
PTRA1590 
PTRA1600 
PTRA1610 
PTRA1620 
PTRA1630 
PTRA1640 


♦  F  2 


M 


ooo  ooo  o  o  o  o  ooo  ooo  o  or* 


f 


i 

i 


IF  (TGI  I)  .EQ.  0.1  Y2( 1 1  *  X6(Ii 

I«IU 

IF  ( I-lBXPl)  270.270,320 
270  IF  (I-ICX-l)  220,220,290 


NEGATIVE  F2  ERROR 


280  Sl-14.0280 
CALL  UNCLE 

NO  SOURCE  IN  ZONE  GREATER  THAN  ICX 

290  IF  (F2.EQ.0.0I  GO  TO  260 
TEMPI l)=H2l I-ll/RRINMUl 
H4I I- 1) *FR EXP (-TEMPI l J-TEMPI ill 
F2*F2*HA( I — 1 I 
GO  TO  260 

300  IF  IF2.EQ.C.0I  GO  TO  310 
TEMPI  1 )*H2I I-ll/RRINMUl 
H4I !-ll*FREXP(-TENP( 11-TEMPI  1) I 
F2*F2*H4I  l~D 
310  5UMX3II)  *  F2 
I»  !♦! 

IF  (I-ICY)  300,300,210 

00  NEGATIVE  ANGLES  SECONO 

320  I  *  1 8XP1 

IF  (IBX-IMJ  370,330.360 
330  IF  CGL!  480,520,340 

GL  »  1/2  HEANS  BLACKBOOY  CONDITION  SET  AT  INP1 
GL  «  POSITIVE  INTEGER  MEANS  INTENSITIES  FROM  QINT1 
GL  *  0  MEANS  VACUUM  AT  IMP1 

GL  NEGATIVE  MEANS  REFLECT IVE  CONOITION  AT  IMP! 

340  IF  (GL.NE.O.SI  GO  TO  350 
F2  »  X6IIMP1I 
GO  TO  480 

350  I QNT  *  Jj  ♦  l  ♦  (NY  ♦  II  *  I IHNU  -  11 
F2  *  ClNTl(lQNT)  /  68.5  *  OHNU 
GC  TO  480 

ERROR  IF  INDEX  EXCEEDS  NORMAL  RANGE 

360  Sl*14.0360 
CALL  UNCLE 

DIFFUSION  BOUNOARY  CONDITION  AT  IBXPl 

370  IF  I  TGI II  .EQ.  0.)  Y2III  -  X6III 
TEMPI  1) =H2( 1 1 /RRINHUI 
IF  I  TEMPI  1 l-i.E-2!  380,380,390 
380  F2X  «  Y2IUU  ♦  TGII  +  1I  ♦  RRINMUl 

F2  *  I(Y2in  ♦  Y2IHIII  *  0.5  ♦  X6I  I)  -  F2X  -  F2X) 
GO  TO  430 

390  H4I  I)*FREXP( -TEMPI  1 1 ) 


PTRA1650 
PTRA1660 
PTRA1670 
PTRA1680 
PTRA1690 
PTRA1700 
PTRA1710 
PTRA1720 
PTRA1730 
PTRA1740 
PTRA1750 
PTRA1760 
PTRA1770 
PTRA1780 
PTRA1790 
PTRA1800 
PTRA1810 
PTRA1B20 
PTRA1830 
PTRA1840 
PTRA1850 
PTRA1860 
PTRA18  TO 
PTRA1880 
PTRA1890 
PTRA1900 
PTRA1910 
PTRA1920 
PTRA1930 
PTRA1940 
PTRA1950 
TABLE  AT  IMP!  PTRA1960 
PTRA1970 
PTRA1980 
PTRA1990 
PTRA2000 
PTRA2010 
PTRA2020 
PTRA2030 
PTRA2040 
PTRA2050 
PTRA2060 
PTRA2070 
PTRA2080 
PTRA2090 
PTRA2100 
PTRA2U0 
PTRA2120 
PTRA2130 
PTRA2140 
PTRA2150 
PTRA2160 
*  TEMPIil  ♦  F2XPTRA21 70 
PTRA2180 
PTRA2190 
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F2=Y2I  I  )*X8(l)  +  (RR(NMU)*TG(I*iI-X8(m*H4<  I) 

GO  TO  430 

399  IF  (TGI  1  +  1 )  .EQ.  0.)  Y2II+1I  »  X6II) 

REGULAR  INTEGRATION  STEP  FOR  F2,  NEGATIVE  HU 

400  IF  I TG(  I )  . EQ<  0.)  Y2( I )  =  X6III 

TEMPI  1) *H2( I I/RRINMUI 
IF  ( TEMPI ll-I.E-2)  410,410,420 

410  F2  =  ((Y2(I»  ♦  Y2I  1  +  1) )  *  0.5  ♦  X6III  -  F2  -  F2I  *  TEMPIil  +  F2 
GO  TO  430 

420  F2=Y2m*X8(I)  +  l(F2-Y2(Ul)-X8I  1*1 1 )*H4I I I+X81 1  +  l I-X8I I )) *H4< I ) 
430  IF  (F2.LT.0.I  GO  TO  460 
440  SUMX4III-F2 

FORM  CONTRIBUTION  TO  X2 

X2I I ) =X2(  I ) -I F2-SUMX3I l ) l*RR( NGSI 

RHOIU-RHOI  I )  ♦(  F2+SUMX3I  I )  )*RR  I NGSI/RRI NMU) 

IF  I  TGI  I )  .EQ.  O.I  Y2I II  -  X6II-1I 
t*t-l 

IF  ( I-IAX)  530,450,450 
450  IF  II-ICY)  500,400,400 
460  Sl=14.0460 
CALL  UNCLE 

NO  SOURCE  IN  ZONE  CESS  THAN  ICY 

470  IF  (F2.EQ.0.0)  GO  TO  480 
TEMPI  1 ) =H2 ( I l/RRINMU) 

H4I  l)=FREXPI-TEMP< 11-TEMPI  1)) 

F2=F2*H4U) 

480  SUMX4I I )=F2 

490  X2 1 1 ) =X2I I »-I F2-SUMX3I I ) I*RR( NGSI 
1  =  1-1 

IF  I I-l-ICX!  399,470,470 

NO  SOURCE  IN  ZONE  LESS  THAN  ICY 

500  IF  (F2.EQ.0.0)  GO  TO  510 
TEMPI  II =H2( I )/RR(NMUI 
H4I 1 1 =FREXP( -TEMPI  1 ) -TEMPI  1 1 1 
F2=F2*H4I U 
510  SUMX4I  I  I=F2 

X2I  IJ  =  X2I  I  )-(F2-SUMX3(  I  HERRINGS) 

1=1-1 

IF  (I-IAX)  530,500,500 
520  F2*0.0 

GC  TO  480 
530  CONTINUE 

IF  I TR08G  .EQ.  0.0  .OR.  TRDBG  .GE.  SOL  10(18)1  GO  TO  539 
C  STORE  INTENSITIES  FOR  DEBUG  PRINT 
XX  *  -0.5 


11 


PTRA2200 
PTRA2210 
PTR',2220 
PTR.  2230 
PTRA2240 
PTRA2250 
PTRA2260 
PTRA2270 
PTRA2280 
PTRA2290 
PTRA2300 
PTRA2310 
PTRA2320 
PTRA2330 
PTRA2340 
PTRA2350 
PTRA2360 
PTRA2370 
PTRA2380 
PTRA2390 
PTRA2400 
PTRA2410 
PTRA2420 
PTRA2430 
PTRA2440 
PTRA2450 
PTRA2460 
PTRA2470 
PTRA2480 
PTRA2490 
PTRA2500 
PTRA2510 
PTRA2520 
PTRA2530 
PTRA2540 
PTRA2550 
PTRA2560 
FTRA2570 
PTRA2580 
PTRA2590 
PTRA2600 
PTRA2610 
PTRA2620 
PTRA2630 
PTRA2640 
PTRA2650 
PTRA2660 
PTRA2670 
PTRA2680 
PTRA2690 
PTRA2700 
P T  RA27 10 
p  TRA2720 
P • R A2  T  30 
Pi KAi  1 40 


jjj  *  JJ  +  1  PTRA2750 

WRITE  (3)  XX,  IAX,  JJJ,  IBXP1,  S0LI0(18),  TH,  RR(NMU)  PTRA2760 


WRITE  (31  ( SUMX31 I ) ,  SUMX4(I», 

XX  =  -2. 

WRITE  (3)  (XX,  I  *  i,  7) 
BACKSPACE  3 

539  DHNU-HNUP-HNU 

540  JJ  *  JJ  +  1 
NMU  *  NHU  ♦  1 
NGS  *  NGS  ♦  1 

IF  (JJ-NYJ  140,140,550 
550  00  560  IMAX.IBXPl 
560  X2(  I)  »=  X2(  I  I*  2.052E12 
RETURN 
END 


IAX,  IBXP1)  PTRA2770 

PTRA2780 

PTRA2790 

PTRA2B00 

PTRA2810 

PTRA2820 

PTRA2830 

PTRA2840 

PTRA2850 

PTRA2860 

PTRA2870 

PTRA2880 

PTRA2890 
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(RAND  Library) 
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313-62,  Project  8814) 

ARMY  ACTIVITIES 
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